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ABSTRACT 

This  thesis  describes  the  study  and  development  of  a 
workable  inverse  scattering  method  for  imaging  and  identifi- 
cation of  radar  targets.   The  space-time  integral  approach 
is  used  for  iterative  target  shape  reconstruction. 

Following  an  overview  of  transient  electromagnetics,  the 
integral  equation  is  applied  for  thin-wire  transient  response 
computation. 

The  analytical  time  domain  integral  equation  is  derived 
and  solved  numerically  for  general  conducting  bodies  of 
revolution.   Finally,  the  algorithm  for  an  inverse  scattering 
computer  solution  is  derived  and  tested  under  simulation  of 
physical  environments. 
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I.   INTRODUCTION 

A.   OVERVIEW 

The  classification  of  radar  targets,  has  become  increasing- 
ly important  in  the  past  two  decades.   In  military  applications, 
the  need  to  distinguish  between  friendly,  enemy,  and  neutral 
targets,  is  vital. 

In  most  radar  applications,  the  only  properties  that  are 
measured  are  range,  angle,  velocity  and  elevation.   The  methods 
that  are  used  by  radars  to  provide  some  classification  of  the 
target  generally  involve  the  examination  of  the  echo  signal. 
Some  of  these  techniques  can  be  summarized  briefly  as  follows: 

1.  High  range  resolution  radar,  which  uses  a  very  short 
pulse  that  can  provide  sufficient  range  resolution  to  obtain 
a  rough  profile  of  the  illuminated  part  of  the  target  shape. 

2.  Cross  section  fluctuations  analysis  which  can  provide 
information  on  the  size  of  the  target. 

3.  Synthetic  aperture  radar  that  gives  a  high  angular 
resolution.  This  method  is  typically  employed  in  ground 
mapping . 

4.  Analysis  of  the  backscattered  polarization  from  a  target 
can  provide  a  means  for  discrimination  between  types  of 
targets . 

5.  Use  of  nonlinear  effects  due  to  junction  of  metals. 
The  backscattered  spectrum  is  analyzed  and  checked  for 


harmonics  appearing  in  the  received  signal,  that  are  typical 
for  certain  type  of  targets. 

6.   Doppler  modulation  by  propellers,  the  rotary  wing  in 
helicopters  or  even  hull  vibrations  in  armored  vehicles. 

None  of  the  methods  cited  can  be  used  to  obtain  all  the 
information  on  the  target  shape  and  size,  they  only  can  be 
used  to  discriminate  to  a  certain  extent  between  different 
types  of  targets.   Other  methods  currently  used  involve 
passive  reception  of  radiation  by  targets. 

Elint  receivers  or  IR  detectors  can  be  used  to  monitor 
radiation  from  targets,  and,  based  on  previous  intelligence, 
or  signature  catalog,  associate  a  type  of  radiation  to  a  type 
of  target. 

The  inverse  scattering  technique , that  is  the  subject  of 
study  in  this  thesis,  uses  the  basic  properties  of  the  transi- 
ent response  to  an  impulse  excitation  (impulse  response)  in 
order  to  obtain  information  on  the  size,  shape  and  also  material 
composition  of  the  target. 

In  principle,  the  smoothed  impulse  response  can  be  ob- 
tained by  measuring  the  backscattered  fields  at  all  frequen- 
cies up  to  some  upper  limit.   In  most  cases  this  would  be 
impractical  so,  instead,  the  target  is  illuminated  by  a  very 
short  pulse  (with  a  very  wide  band  frequency  spectrum  ranging 
from  the  low  frequencies  in  the  Rayleigh  region  up  to  higher 
resonance  region  frequencies) . 
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The  impulse  response  contains  all  the  information  needed 
to  characterize  the  target.   The  higher  frequencies  give  the 
information  about  the  structure.   This  is  related  to  the 
optical  theory  of  diffraction  which  is  a  good  approximation, 
under  certain  conditions,  for  shape  prediction  up  to  the 
shadow  boundary.   Physcial  optics,  however,  fails  to  predict 
the  phenomena  appearing  in  the  shadow  region,  such  as  the 
creeping  wave  which  is  due  to  lower  frequency  waves  travelling 

around  the  body . 

A  different  approach  that  has  been  used  to  retrieve  tar- 
get configuration  information  is  related  to  the  singularity 
expansion  method  (S.E.M.).   This  method  seeks  to  classify 
scatterers  through  the  complex  residues  and  poles  of  their  time 
domain  signatures.   Briefly,  in  this  method  the  transient  E.M. 
response  of  a  conducting  body  is  characterized  as  a  series  of 
complex  exponentials  (or  damped  sinusoids):  [1] ,  [2] 

N      st 

f(t)    =     z»  .- 

m=l 

where  R  is  the  amplitude  (residue)  of  each  mode  of  complex 
m 

frequency  s  .   The  true  function  has  a  complex  frequency 
representation  of  the  form: 

N   R 
F(s)   =    I 


m=l     m 


The  poles  s   depend  only  on  the  geometry  (independently 
r      m 

of  the  orientation) ,  and  the  residue  R   depends  on  both 
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excitation  and  geometry  [3].   Thus  once  the  poles  and  residues 
are  known,  they  provide  a  means  for  target  classification. 
The  problem  is  to  find  the  poles  and  residues. 

A  method  is  suggested  using  Prony's  algorithm  in  [2]  for 
extraction  of  poles  and  residues  from  time-domain  signatures. 
The  method  has  been  tried  successfully  on  wire  geometries  and 
simple  shapes. 

Another  important  application  of  the  impulse  response  of 
a  target  is  the  study  of  its  radar  cross  section  (RCS) ,  since 
the  Fourier  transform  of  the  impulse  response,  (the  frequency 
response)  is  directly  related  to  the  RCS  by  the  following 
relationship: 


2  i-/  »  i2 
a   =   it  c    F  (co) 


where 


F(oj)   =     /  F  (t)  eZ,a)t  dt 


I 

—  CO 


and 


F  (t)  is  the  impulse  response  of  the  target. 

Finally  the  time-domain  analysis  of  the  impulse  response  of 
a  body  is  of  great  importance  in  the  study  of  the  electro- 
magnetic pulse  (E.M.P.)  response  of  targets. 
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B.  PRACTICAL  CONSIDERATIONS 

The  practical  study  and  measurement  of  the  impulse  re- 
sponse involves  the  building  of  a  time-domain  scattering  range. 
The  availability  of  generators  having  very  short  pulse-widths, 
fast  rise  time,  and  high  voltage  outputs,  together  with  fast 
sampling  oscilloscopes  and  small  computers  makes  the  task 
feasible . 

For  the  solution  of  the  inverse  scattering  problem  the 
time  domain  scattering  range  developed  at  the  Naval  Post- 
graduate School  [4] ,  was  employed  to  measure  the  impulse 
response  and  demonstrate  the  validity  of  the  theories  developed 
in  this  thesis.   The  transient  response  data  measured  by  using 
pulse  generators  with  very  fast  rise  times  can  validate  the 
calculation  of  the  time  domain  integral  equation  solution. 
More  important,  it  provides  a  vehicle  for  testing  various  tar- 
get imagery  and  I.D.  schemes,  as  well  as  investigating  broad- 
band RCS  reduction. 

C.  SURVEY  OF  THE  LITERATURE 

One  of  the  first  efforts  to  determine  physical  properties 
of  electromagnetic  scatterers  from  time-domain  analysis  was 
due  to  Kennaugh  and  Cosgriff  [5] ,  who  employed  the  physical 
optics  approximation  to  calculate  the  approximate  backscattered 
impulse  response  of  a  flat  spheroid. 

A  summary  of  this  work  was  reported  by  Kennaugh  and  Moffatt 
[6] ,  who  extended  the  physical  optics  approximation  to  transient 
response  calculations.   At  this  stage  the  computation  of  the 
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impulse  response  was  approximated  by  using  the  physical  optics 
approach  or  the  inverse  transform  of  the  frequency  domain  solu- 
tion.  Then,  in  1968,  Bennett  derived  in  his  Ph.D.  thesis  [7] 
an  approach  for  solving  the  time  domain  integral  equation 
directly  by  using  a  time -stepping  method.   Bennett  used  a 
form  of  the  equation  called  the  magnetic  field  integral  equation 
(or  M.F.I.E.)  for  two  and  three-dimensional  surfaces. 

Sayre  and  Harrington  [8]  used  the  same  time-stepping  tech- 
nique to  derive  the  transient  response  of  a  wire,  but  they 
employed  another  type  of  equation  called  the  electric  field 
integral  equation  (EFIE) . 

Additional  work  in  time-domain  electromagnetics  has  been 
reported  in  several  more  recent  papers  as,  for  example,  the 
report  by  Miller,  Poggio  and  Burke  [9]  on  time-domain  analysis 
of  thin-wires.   Summaries  and  review  of  the  time-domain  electro- 
magnetics are  due  to  Mittra  [2]  and  [10]  and  Miller  [11] . 

A  comprehensive  review  on  the  subject  and  its  application 
has  been  published  by  Bennett  [12],  which  gives  the  reader 
an  introduction  and  initiation  to  time-domain  electromagnetics. 
Literature  dealing  with  the  inverse  scattering  problem  is  due 
to  Young  [13]  and  Shubert,  Young  and  Moffatt  [14]  ,  which  des- 
cribes the  image  generation  of  a  target  from  the  ramp  response 
waveform. 

Moffatt  and  Mains  [15]  describe  a  method  of  detection  and 
discrimination  of  radar  targets  using  the  idea  of  ramp  response 
but  with  multiple  harmonic  radar  frequencies.   Finally,  a  very 
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recent  report  by  Bennett  [16]  gives  the  general  approach  to 
the  time  domain  inverse  scattering  solution.   The  method  sug- 
gested in  that  report  is  being  used  in  this  thesis  to  approach 
the  problem  of  transient  target  imaging.   This  list  of  efforts 
is  by  no  means  exhaustive,  as  many  other  authors  have  published 
papers  dealing  with  transient  electromagnetics  and  target 
imaging. 

D.   THESIS  OBJECTIVES 

,The  objective  of  this  effort  is  to  develop  the  analytical 
and  numerical  solution  of  time-domain  inverse  scattering 
problem.   It  follows  the  development  of  the  time-domain 
scattering  range  at  the  NPS  [4]. 

The  main  effort  has  been  conducted  toward  the  formulation 
and  numerical  solution  of  the  time-domain  integral  equations 
for  conducting  bodies  of  revolution. 

A  computer  program  is  developed  to  solve  for  direct  and 
inverse  scattering.   The  algorithm  was  used  for  experimental 
display  of  the  shape  of  the  target  under  test. 
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II.   GENERAL  SOLUTION  OF  THE  TIME 
DOMAIN  INTEGRAL  EQUATIONS 

This  chapter  is  a  review  of  the  development  of  time  domain 
integral  equations,  their  solutions  and  the  associated  numeri- 
cal considerations.   Detailed  development  of  the  equations 
can  be  found  in  [2,10,11]. 

A.   THE  INTEGRAL  EQUATIONS  (E.F.I.E.  and  M.F.I.E.) 

The  basis  for  derivation  of  the  equations  are  the  general 
time  dependent  forms  of  Maxwell's  equations: 


V  xE(r,t)   =   -  ft  yo  **<r,t) 


V  xH(r,t)   =   |r-  £   E(r,t)  +  J(r,t) 

ot   O 


V  •  E(r,t)  =   p(r,t)/eo 


V  •  HQ(r,t)   =   0 


Starting  from  these  equations,  and  using  the  expressions  for 
the  magnetic  vector  potential,  A,  and  the  electric  scalar 
potential,  $ ,  which  are  related  to  the  fields  by 


H(r,t)   =   —  7  xA(r,t) 
^o 


E(r,t)   =   -  7  (j>(r,t)  -  |£(r,t) 


16 


one  can  obtain  the  expression  of  the  magnetic  field  integral 
equation  (M.F.I.E.)  for  a  perfect  conductor  [2] 


Je(r,t)   =   2n  xHinc(r,t) 


+  2T^X/  Cc-|T  +  I]5s(^^  *4dS'      (2.D 

D  R 

and  the  electric  field  integral  equation  (EFIR) 

->  II 

nxE         (r,t)       =      _x/[__  Js(r',x) 

s 
.  £i£^il(i  +  iLJids-  (2.2) 

£    RZ         R         C    3T 
O 

In  equations  (2.1)  and  (2.2), 


n   is  the  unit  vector  normal  to  the  conductor  at 
the  observation  point  r. 

-y  .     ->■ 

J  (r',x)   is  the  current  density  at  the  source  point  r1 
at  a  retarded  time  x . 

o(r',T)   is  the  charge  density  at  source  point  at  the 
retarded  time  x . 


Here  x  is  given  by 

x   =   t  -  R/c 


and 


R   is  the  distance  from  the  observation  point  r  to 
the  source  point  r',  i.e., 


R   =   r  -  r ' 
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->inc  ■*■   x      -^inc  -*■ 

E    (r,t)  and  H    (r,t)  are  the  incident  electric  and  magnetic 

fields  at  the  observation  point.   Because  of  numerical  con- 
straints, the  equation  most  suitable  for  closed  surfaces  is 
the  MFIE  equation  (2.1)  while  for  thin  wires  and  open  struc- 
tures the  EFIE  equation  (2.2)  is  used. 

Equation  (2.1)  states  that  the  current  induced  on  a  point 
r  of  a  body  at  a  time  t  is  due  to  the  incident  field  H    at 
that  point  and  the  sum  of  scattered  fields  there  from  earlier, 
more  distant  locations  of  the  body.   The  interactions  between 
the  samples  of  the  body  are  displaced  in  time  by  an  amount 
equal  to  the  transit  time  of  the  fields  between  them. 

In  the  integral,  the  observation  point  (R  =  0)  is  excluded 
since  this  "self-term"  is  included  in  the  incident  field  por- 
tion of  the  equation.   This  determines  the  general  solution 
for  the  MFIE  equation  which  is  solved  directly  by  marching 
on  in  time.   The  solution  is  started  at  a  certain  time  (t  =  0) , 
and  due  to  causality  the  integral  part  of  the  equation  will 
be  zero  since  there  are  no  past  currents.   Then,  as  t  increases, 
the  integral  is  solved  with  previous  currents  already  computed. 
This  technique  is  a  time-stepping  procedure  for  solving  initial 
value  problems. 

For  the  EFIE,  equation  (2.2),  the  method  is  slightly  differ- 
ent but  uses  the  same  time  stepping  procedure.   The  solution 
of  this  equation  is  applied  to  the  problem  of  transient 
scattering  by  a  thin  wire,  in  the  next  chapter. 

The  solution  of  equations  (2.1)  and  (2.2)  gives  the  cur- 
rents induced  on  the  body  by  the  incident  field.   Once  the 
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currents  are  known,  it  is  a  straightforward  step  to  compute 
the  scattered  fields.   For  the  far-field  case  the  field  is 
given  by : 


$scatt  ,-> 
n 


(?o't}    ■    4iF-5-   /  I?  5(r,'T)  x*Rds'  (2-3> 

O   S 

where  aR  is  the  unit  vector  from  the  source  r '  to  the  far 

observation  point  and  r   is  the  separation  distance.   J  (r',x) 

is  the  current  at  the  source  point  at  a  previous  (retarded) 

time  x  =  t  -  r  /c. 
o 

B.   NUMERICAL  SOLUTION  AND  CONSIDERATIONS 

The  solution  of  the  integral  equations  implies  finding 
the  currents  induced  on  the  body  once  illuminated  by  an  inci- 
dent field.   Since  the  solution  is  to  be  carried  out  numerically, 
the  first  step  is  to  represent  the  currents  by  their  space  and 
time  samples  J. .(s,t)  in  a  discrete  manner.   J. .(s,t)  is  the 
current  at  the  i    segment  AS  (or  patch)  at  the  j    time 
increment  AT.   The  total  current  will  be  approximated  in  terms 
of  all  samples  at  all  times  and  spatial  segments  via  the  basis 
function  expansion: 


NS   NT 


J(r,t)   =    I         I      J.  .  (s,t)  U(t  .)  V(s.  : 
i=l  1-1   ^  1 


(2.4) 


where  N„  is  the  number  of  space  samples  on  the  body  while  NT 
is  the  total  number  of  time  increments,  and  the  pulse  basis 
functions  are: 
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U(tj)   = 


v(si)  = 


1  at  the  center  of  the  jth  time  interval 

0  elsewhere 

1  at  the  center  of  the  ith  space  patch 

0  elsewhere. 


The  space  time  sampling  is  shown  in  figure  2.1.   The 
sampling  of  the  currents  already  suggests  that  the  body  is 
divided  into  segments  (or  patches)  AS  and  the  time  is  sampled 
at  AT  intervals,  while  the  total  observation  time  goes  up  to 
NT  x  AT. 

It  is  very  important  at  this  point,  before  attempting  to 
solve  this  equation,  to  define  the  basic  criteria  for  proper 
selection  of  the  space  and  time  sampling  as  well  as  their 
relationship  to  the  shape  of  the  exciting  pulse  and  the  body 
size. 

The  incident  field  is  chosen  here  to  be  represented  by  a 
gaussian  impulse  due  to  its  rapidly  decreasing  features  and 
also  because  most  of  the  short  pulse  generators  will  produce 
such  a  waveform. 

The  gaussian  impulse  is  given  by 

g(t;   =   e 

Its  amplitude  falls  to  1/10  of  its  maximum  value  at 

t     -  1.5/A.   Its  frequency  spectrum  is  given  by  [2] 
max 

2  2   2 

G(f)   =   F{g(t)  )   =   e-TT  r  /A  . 
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computation  point 


Figure    2.1.      Space-time    sampling   diagram 
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This  response  reduces  to  1/10  of  its  maximum  value  at 
fmax  ~  0,5A-   The  time  sampling  AT  must  exceed  the  Shannon 
rate  to  avoid  aliasing  problems  if  FFT's  are  employed  later, 
i.e.  , 


AT     <_     jzr —  _   1 
nun       max     A 


Usually  AT  will  be  taken  smaller, 


AT 


5fmax   =   275A  (2'5) 


The  spatial  sampling  AS  must  be  chosen  such  that 

AS   >   AT  xc  (2.6) 

otherwise  there  will  be  coupling  between  the  currents  of 
adjacent  patches.   Violation  of  (2.6)  will  also  mean  that  the 
stepping  procedure  for  solving  the  equations,  which  assumes 
that  during  one  time  interval  the  currents  at  one  patch  are 
not  influenced  by  adjacent  currents,  will  no  longer  be  satis- 
fied, thus  producing  an  incorrect  solution. 

Another  consideration  is  that  the  space  sample  points  must 
be  close  enough  to  resolve  adequately  the  spatial  variation 
of  the  incident  field  over  the  body.   From  Eq.  (2.5)  and  (2.6) 
we  have : 


AS   > 


5f 
max 
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A  . 

as  >   mi-n 


This  latter  relationship  imposes  some  minimum  ratio  of  the 

size  L  of  the  body  to  the  minimum  wavelength  (or  equivalently 

the  minimum  pulse  width) .   The  ratio  L/A  .   determines  the 

min 

resolution  and  the  usefulness  of  the  time  domain  analysis, 
and  it  must  be  kept  high  enough  to  maintain  adequate  resolution 
If  the  body  is  to  be  represented  with  good  fidelity  with  a 
certain  minimum  number  of  samples  N  we  will  need 

Ne  A  . 
T   ,    S   mm 
L   >_  g 

This  means  that  if  A    increases  (i.e.,  wider  pulse)  L  must 

mm  v  ' 

increase  or  equivalently  for  a  small  size  target  A  .   must 

^         J  r     min 

decrease  (i.e.,  shorter  pulse). 

Having  presented  the  equations,  their  general  solution  and 
the  requirements  for  the  various  parameters  involved  we  can 
move  on  to  the  implementation. 

C.   APPLICATION  TO  THE  TIME -DOMAIN  SCATTERING  BY  A  THIN-WIRE 

The  first  application  that  was  studied  was  the  computation 
of  the  impulse  response  of  a  straight  thin  wire.   This  example 
was  chosen  since  it  gives  a  good  understanding  of  the  physical 
processes  and  the  numerical  techniques  used  in  the  solution 
of  time  domain  integral  equations.   The  basis  for  this  compu- 
tation is  given  by  Sayre  and  Harrington  in  [8]. 

The  problem  is  formulated  as  follows:   a  thin-wire  of 
length  L  and  radius  a  aligned  with  the  z  axis  is  illuminated 
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by  an  impulse  field  of  the  form 

-A2(t-t    )2 
Einc(t)   =   e       max 


The  basic  formulas  are  in  terms  of  the  retarded  potentials : 
The  boundary  condition  : 

Cc<^>  -'t^  +  H 

where  the  retarded  magnetic  potential  is  given  by  : 

a  (t  Z)   =  JL     (   I(z' /T)  dzi 
Az(Z,Z}  4tt  r  .  1  R    dz 

Axis 

and  the  retarded  electric  potential  is 

<0   =   — ±—  \    q(2'>T)  dz' 

o  Axis 

where  the  continuity  equation  between  charge  and  current  is 


91        3g 

8z        3t 


In  these  equations  I(z',t)  and  q(z',T>  are  the  current  and 
charge  at  source  point  z'  at  time  x  =  t  -  R/c ,  while  R  is 
the  distance  between  the  source  and  observation  points .   The 
geometry  is  described  in  figure  2.2. 

Using  the  numerical  formulation  expressed  in  [8],  and 
applying  the  time  stepping  procedure,  a  computer  program  was 
written  (Appendix  D)  to  compute  the  scattered  field  from  a 
thin  wire  of  length  1  meter  and  radius  0.8  mm.   The  wire  was 
divided  into  40  segments  and  the  sampling  time  AT  =  AZ/c  =  L/40c 
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Figure  2.2.   Geometry  of  the  thin-wire  scatterer 
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The  incident  field  was  taken  to  be  similar  to  the  pulse 
transmitted  by  the  generator  of  the  NPS  transient  scattering 
range,  with  linear  polarization  parallel  to  the  wire  and 
incidence  normally  on  it. 

2         2 
-A^(t-t    )' 
„inc  max 

E      =   e 

A  =   6 »10   sec 


t     =   18-AT  (about  half  transit  time  along 
max  the  wire) 


The  program  computes  the  currents  at  each  spatial  segment  with 
up  to  400  time  samples.   This  corresponds  to  10  transit  times 
across  the  wire.   The  current  at  the  center  of  the  wire  is 
given  in  figure  2.3.   The  backscattered  field  is  given  in 
figure  2.4. 

No  significant  change  was  obtained  for  20  segments  and 
200  time  increments.   Comparison  of  this  computation  with  lab 
measurement  was  done  by  Hammond  [4]  and  very  good  agreement 
was  found. 

Similarly,  comparison  of  the  computed  results  found  in 
[11]  with  a  ratio  of  L/a  =  0.00667  and  A  =  3.25-10   gave 
exact  agreement.   The  backscattered  field  for  this  case  is 
shown  in  figure  2.5. 
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Figure    2.3 


Current   induced    at    the    center   of    a 
thin-wire   by   a   gaussian    incident  pulse 
(A  =    6-108    and    a/L   =    0.0008) 
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Figure  2.4.   Backscattered  field  from  a  thin-wire 
(A  =  6-108,  a/L  =  0.0008) 
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Figure    2.5.       Back scattered    field    from   a    thin-wire 
(A   =    3.25-109    and    a/L    =    0.0067) 
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III.   TIME  DOMAIN  RESPONSE  OF  PERFECTLY 
CONDUCTING  BODY  OF  REVOLUTION 

A.  INTRODUCTION 

This  section  considers  the  solution  of  the  time  domain 
response  for  a  body  of  revolution,  when  the  exciting  wave  is 
propagated  along  the  axis  of  symmetry  of  the  body. 

The  use  of  the  axial  symmetry  of  the  body  will  result  in 
considerable  simplification  of  the  computation  of  both  the 
induced  currents  on  the  body  surface  and  the  field  scattered 
from  the  body.   The  algorithm  developed  in  this  section  will 
be  used  later  for  the  solution  of  the  inverse  scattering 
problem. 

B.  THE  MAGNETIC  FIELD  INTEGRAL  EQUATION 

In  [7]  Bennett  derives  the  integro-dif ferential  equation 
for  three  dimensional  scatterers: 


-*■-*■  ~    +i  ->■ 

J(r,t)   =   2a  xH  (r,t) 
n 


+  ^^n*1^^!?'^'^*!1^'     (3-X) 


where 


x  =   t  -  R/C  (retarded  time) 

-> 

r  is  the  observation  point 

r*  is  the  source  (or  integration)  point 
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R  =   |r  -  r1  I 

A  is  the  surface  of  the  scatterer 

"*"i  ■*■ 

H  (r,t)   is  the  incident  magnetic  field  at  the  observer's 

point. 

The  geometry  for  the  body  of  revolution  is  given  in  figure 
3.1. 

For  the  actual  problem  the  incident  magnetic  field  is 
propagating  along  the  z-axis  (axis  of  rotation  of  the  body) , 
and  is  directed  towards  the  positive  x-axis.   The  body  lies 
completely  to  the  right  of  the  positive  z-axis,  and  is  defined 
by  the  continuous  contour  function,  p(z),  which  is  the  radius 
of  the  body  at  a  point  z. 

The  prime  coordinates  denote  the  source  or  integration 
point,  while  the  unprimed  denote  the  observation  point  for 
which  the  current  is  desired.   At  each  point  on  the  body,  there 


/\     /\ 


is  defined  the  triad  of  orthogonal  unit  vectors,  (a  ,a  , a  ) 
where : 

a  :   unit  vector  normal  to  the  surface 
n 

a , :   unit  vector  tangent  to  the  surface  and  directed 
'    towards  the  positive  <f>  direction 

a  :   unit  vector  tangent  to  the  surface  and  directed 
along  the  curve  representing  the  longitudinal 
contour  of  the  body. 


For  the  following  derivation  we  will  define  the  angle  a  (or 

a')  as  that  sustained  by  a   (or  a')  with  the  z  axis,  a  (or  a') 

s       s 

being  positive  when  a   (or  a')  points  outward  from  the  z  axis 

S  S 
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and  negative  if  it  points  towards  it.   The  angle  tj>  (or  <J)  ' ) 
is  the  angle  defined  by  p  (or  p1)  with  the  (x,y)  plane. 
Equation  (3.1)  states  that  the  current  density  at  the  obser- 
vation  point  r  at  time  t  is  due  to: 

1.  The  current  induced  by  the  incident  field  at  (r,t) . 

This  part  represents  the  physical  optics  approximation, 

2a  xUx(rrt). 
n 

2.  The  current  due  to  the  flow  of  currents  at  all  other 
points  of  the  body  at  a  retardated  time  t  =  t  -  R/C,  where 
R  is  the  distance  from  the  integration  point  to  the  oberva- 
tion  point. 

According  to  the  geometry  described  in  figure  3.1  we  can  write 

R   =   (p2  +p'2  -  2pp'  cos  (<f>'  -(j))  +  (z'  -z)2)1/2 

The  current  on  the  body  will  flow  only  in  the  s  and  <p    direc- 
tions and  can  be  expressed  (for  both  primed  and  unprimed 
coordinates) ,  as 

J(r,t)   =  Js  (p,z,t,<j>)  -as  +  J  (p,z,t,c|>) -a, 

where  J  (p,z,t,cj))  and  J,  (p ,  z,  t,<j>)  are  the  scalar  values  of 
the  current  components  and  are  functions  of  the  position  on 
the  scatterer.   We  will  denote  them  simply  by  J   and  J,  so 
that 

J(r,t)   =   J  ao    +   J,  a,  . 
s   s     (p   <p 


32 


Figure  3.1 


Geometry  of  the  scattering  by 
a  body  of  revolution 
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Eq.  (3.1)  can  be  written  in  terms  of  its  components  in  the 
s  and  $  direction,  i.e., 


J(r,t)   =   [2an  xH1(r,t)].  +i  ji[i  +  i  |-]  [a  x  J-  x  R]„  dA' 

s  n         s   ztt  -.J  ^z  R   c  dx   n        a 

63      A  R  b 

(3.2) 

J.(r,t)   =   [2an  xH1(r,t)]^  +i  /  *  [i  +  i  |_]  [£  x  j  '  x  R^dA' 

(3.3) 

In  order  to  obtain  the  scalar  expressions  of  the  two  last 
equations  we  have  to  perform  the  vector  multiplications: 
This  is  done  in  Appendix  A  and  as  a  result,  equations  (3.2) 
and  (3.3)  become: 

Js  =  -2^(8^)8111  ♦  +A-  /i[i  |7  +  |][J'.F1-J'V.sing]dA- 

A  R  r 

(3.4) 

1     flrl   3    .1. 


J.   =   -2H1(z,t)cos  <J>  sin  a+^p-  f-=k-[=r  i—  +  =-][J'  U  sin  3  +J'F«]  dA1 

<p  2tt  a  r^  C  3t   R    s  $  2 

(3.5) 
where : 

3  =  <j>'  -  <t> 

F,  =  (z'-z)cosBsinc^  +  (p  -  p  *  cos  3  )  cos  a  ' 

V  =  z1  -  z 

F„  =  (p  cos  3  -  p' )  cos  a  +  sin  a  cos  3  (  z  '  -  z) 

U  =  (z  '  -  z)  sin  a  •  sin  a'  +p  cos  a  sin  a'  -p'sina-cos  a' 
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,1     . 


We    note    that    the    incident   field   H      is   dependent   only   on    z   and 
t    (independnet   of    <J>    due    to    symmetry)  .       We    also   note    that: 


dA' 


p'    d<p'    ds' 


p'    dB   ds' 


The  integrals  in  Eq.  (3.4)  and  (3.5)  can  then  be  divided  into 
two  parts:   integration  over  3  from  0  to  2tt,  and  integration 
over  ds'  from  0  to  S  (S  being  the  contour  curve  as  shown  in 
figure  3.2). 

♦  P  (z) 

S 


Figure  3.2.   Body  contour 

The  next  step  is  to  obtain  further  simplication  of  the  inte- 
gral equations  (3.4)  and  (3.5). 

The  currents  in  these  equations  are  functions  of  their 
position  on  the  body: 


Js(z,p,<|) ,  t) 


J    =   J  (z,p,<J>,t) 

But  we  remark  that  the  currents  due  to  the  incident  field  have 
a  sinusoidal  variation  with  <j> ,  as  follows: 


J    =   H  (z,t)  cos  (f  •  sin  a 
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J        =      H    (z,t)     sin   <j> 

so   that  we   can   assume    that   the  currents  will  have    the    same 
respective    relationships  with  <j) . 

J    (z,p,t,<f>)       =      J    (z,p,t)sin   (()      =      J    -sin    <J>  (3.6) 

o  S  5 

J,  (z,p,t,(j))       =      J,(z,p,t)cos    d>      =      J    -cos    d>  (3.7) 

9  <P  cf> 

This  assumption  is  demonstrated  in  Appendix  B,  using  Maxwell's 
equation. 

After  substitution  of  equations  (3.6)  and  (3.7)  into  the 
integrals  (3.4)  and  (3.5),  the  integrals  over  the  angle  3  can 
be  simplified  as  shown  in  Appendix  C,  so  that  the  scalar 
integral  equations  reduce  to  their  final  expressions: 

•[J'F,  cos  3  +Vsin2BJ']d$  (3.8) 

si  tp 


J    =   -2Hi(z/t)sina+  A.  J    p .  ds  -   J    1[|  +  1|_] 


•[J1  U  sin23  +  F„  cos  3  J ! ] d3  (3.9) 

S  2  <p 

in  which  the  currents  J   and  J,  (J1  and  J!)  are  only  functions 

S        (J)     s        (J)  J 

of  p ,  z,  and  t  (p'  ,z'  ,t)  and  independent  of  the  angle  <j> . 

The  significance  here  is  that  the  surface  distribution  of 
the  unknown  currents  is  reduced  to  a  one-dimensional  variation 
of  current  along  the  generating  arc  of  the  surface  of  revolution 
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This  will  allow  considerable  reduction  in  computer  time  and 
storage  when  carrying  out  the  numerical  solution  of  the 
equations . 

Once  the  currents  along  the  contour  of  the  body  are  known, 
the  currents  at  any  other  point  (at  some  time  t  and  some  z 
position)  will  be  given  by 


Js(p,z,t,<j>)   =   Js(p,z,t)  -sin  <J) 


J  (p,Z,t,4>)    =   J,  (p,Z,t)  'COS  c)> 


C.   FAR -ZONE  SCATTERED  FIELD 

In  most  cases  we  are  interested  in  the  far  scattered 
field.   At  a  large  distance  from  the  body  the  scattered  field 
is  given  by: 


HS(?'t}   =   4^TC   /fr^r',!)  x£R]dA'  (3.8) 


o   A 


where 


x   is  the  re tar dated  time  t  -  R'/C 

R'   is  the  path  difference  between  the  integration 
point  to  the  observer  and  the  origin  to  the 
observer. 

-»-  ->- 

J(r',x)   is  the  current  at  the  surface  integration  point 

retarded  in  time,  while  a   is  the  unit  vector 
pointing  from  the  source  point  towards  the  obser- 
vation point.   The  geometry  is  described  in 
figure  3.3. 

r   is  the  distance  from  the  origin  to  the  observation 
o      . 
point. 
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direction  of 
propagation 


A  is   the   projection   of    the   source   point   on   the 

plane    (x,z) 


Figure  3.3.   Geometry  for  the  far-zone 
scattered  field 
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In  the  following,  the  scattered  field  will  be  derived 
for  an  arbitrary  elevation  angle  6.   Due  to  symmetry  with 
<j>  there  is  no  loss  of  generality  if  we  observe  the  scattered 
field  in  the  plane  (x,z)  where  <J)  =  0 . 

According  to  the  geometry: 

R'   =   ~[z'  cos  8  +  p'  cos  <j> '  sin9] 
x   =   t  -  tf/C. 

R'  is  negative  whenever  the  observation  point  is  closer  to  the 
source  than  the  origin,  and  positive  otherwise. 

a„   =   sin  0  a„  +  cos  9  a 
R  x  z 

Performing  the  vector  multiplication  in  equation  (3 . 8)  and 
separating  the  field  into  its  cartesian  components,  gives  the 
f ollowoing  set  of  three  scalar  equations : 


(3.9) 


H      6,t        =      .        n      /   I— [J   sin   4>' sin   a'+J.cosV    cos  8dA' 
x      '  4iTr   C   _ J     3x      s  y  d>  T 

O      A  r 

H    (9  ,  t)       =      -. ^r      f    7 — [  (J  .   -  J     sin  a  ' )  sin   <j>*  cos    (f)1  cos    0 

y  4irr   C    ,;     ox         tb         s  Y  T 

2  O      A  Y 

+   J      sin   (J)'  sin  9  cosa']dA'  (3.10) 

HS(9,t)       =       .    -    _      /   i-[-  sin  9  (JsinVsin   a'  +JAcos2<}) '  ]dA' 
z  4irr   C    „ ;     3x  s  Y  d> 

o      A  Y 


(3.11) 


After  substitution  of  dA'  =  p '  d  d/  ds  '  ,  and  using  the  odd 
properties  of  some  of  the  terms  appearing  in  the  integral  on 
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<j> '  (as  it  has  been  done  previously  in  Appendix  C)  ,  the  three 
scalar  equations  representing  the  components  of  the  scattered 
field  reduce  to  the  following  expressions: 

HS(0,t)   =   cos  9  HS(t)  (3.12) 

X  o 

HS(e,t)   =   -  sin  9  HS(t)  (3.13) 

z  o 

HS(9,t)   =   0  (3.14) 

g 

where  the  common  term  H  (t)  is  given  by  the  following  integral: 

Ho(t)   =  Ithc   J    p'ds'  J   *  fr(Js  sin2*'  sin  a- 

+  J  cos^'W  (3.13) 

Equations  (3.12)  to  (3.15)  represent  the  scattered  field  for 
an  arbitrary  elevation  angle  9 .   The  special  case  of  interest 
is  the  backscattered  field  where  9  =  it  .   For  this  case  there 
is  only  one  component  of  the  scattered  field,  which  is 

H*(t)   =   -  H«(t) 

and  the  retarded  time  is: 

t   =   t  -  R'/C 

where 

R*   =   z'. 
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All  the  terms  appearing  inside  the  integral  in  equation 
(3.15)  are  independent  of  <j> '  so  that  the  integration  on  <j> ' 
can  be  carried  out  analytically.   Doing  this  yields  the 
following  expression  for  the  backscattered  field: 

i     s  a 
Hs(t)   =  ~  IT~C      $      J7[J  S(p'/Z'/T)sin  a'  +J  (p  '  ,  z  '  ,t)  ]  p 'ds  ' 
o   0 

(3.16) 

In  equation  (3.16),  the  currents  J   and  J,  are  those  given 
by  equations  (3.8)  and  (3.9)  respectively,  used  at  time 
t  =  t  -  z*/C  (z1  >_  0)  ,  where  z'/C  is  the  retardation  time  rela- 
tive to  the  origin.   The  backscattered  field  equation  is 
independent  of  <j)  which  means  that  we  can  integrate  the  time 
derivative  of  the  currents,  along  a  ]ine  on  the  contour  of  the 
body. 

D.   NEAR  ZONE  BACKSCATTERED  FIELD 

In  the  Time  Domain  Scattering  Range  developed  at  NFS, 
the  distance  from  the  antenna  to  the  target  under  test  is  not 
long  enough  to  be  considered  for  all  frequencies  of  interest 
as  a  far  region.   Consequently,  it  is  needed  to  develop  a 
formulation  of  the  scattered  field  in  the  near  zone.   For  this 
case  the  geometry  is  described  in  figure  3.4. 

The  general  equation  for  the  scattered  field  is  given  by: 


SS(?,t)    -   ^/l-^l-,  [?<?■, T,  x8]*l  <3.17) 
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direction  of 
propagation 


observation 
point 


Figure  3.4, 


Geometry  of  the  near-zone 
backscattered  field 
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where 


t   =   t  -  R'/C,  and 
dA'   =   p'  df  ds1 

R2  =   (r2  +  z'2)  +  p'2 
o 

After  carrying  out  the  vectorial  multiplication,  dividing  the 
field  vector  into  its  cartesian  components,  again  dropping 
all  the  odd  terms  in  cj> '  and  integrating  analytically  on  <£ '  , 
equation  (3.17)  reduces  to: 

Hx(ro'tJ    =   k/  ^72cR^  +  ?lr][Js(p,cosa,-(ro  +  z,)sin  a' 

-  J.  (r  +z,)]p,ds' 

*   °  (3.18) 

Equation  (3.18)  is  the  expression  of  the  near  backscattered 
field  where: 


((r  +z')2  +  p,2)1/2 


The  currents  J  (p',z',t)  and  J, (p',z',x)  are  given  by  equations 
(3.8)  and  (3.9),  being  used  at  the  retard  time  x. 


t   =   t-  (R'-r  )/C    (referred  to  the  origin) 
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IV.   NUMERICAL  SOLUTION  OF  THE  TIME  DOMAIN 
RESPONSE  OF  THE  BODY  OF  REVOLUTION 

In  the  previous  section  the  analytical  equations  for  the 
currents  and  scattered  field  were  derived.   The  equation  ob- 
tained will  have  to  be  solved  numerically,  and  the  first  step 
will  be  to  divide  the  body  into  patches  of  approximately  equal 
curvilinear  surface  area.   This  is  followed  by  the  setup  of 
the  algorithm  for  the  solution. 

A.   DIVIDING  THE  BODY  INTO  PATCHES 

The  geometry  of  the  body  is  completely  specified  by  the 
function  p(z),  (0  <  z  <_  L) ;  which  represents  the  body  contour 
function  S.   Since  the  integration  is  to  be  carried  over  S, 
the  contour  will  be  divided  into  NS  segments  of  equal  length 
AS.   Consequently,  the  body  contour  function  is  fully  defined 
by  NS+1  points. 

p ^ { z  )       I   =   1   to   NS+1 

The  space  samples  of  each  segment  AS  is  at  the  middle  of  the 
segment  and  is  given  by 


Pi   =   (P£+1  +  P£)/2 


£  =  i  =  1, . . . ,NS    (4.1) 
Zi   =   (Z£+1  +  Z£)/2 
The  angle  a  for  each  sample  point  is  the  slope  of  the  curve 
at  that  point  and  is  approximated  linearly  by 
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P£+l  "  P£ 
a.   =   Arctan  (-Zl± i)  f   £  =  i  =  1,...,NS       (4.2) 

1  Z£+l  ~  Z£ 


Figure  4.1  shows  the  geometry  of  the  contour  sampling.   Each 
segment  AS  defines  a  ring  at  z.  with  radius  p. .   All  the  rings 
will  be  divided  into  patches  of  approximately  equal  surface 
area,  dA.  ,  according  to  the  following  f ormulas j 

dA.   =   DS  •  p.  •  AS.   =   AS2 
l  l      l 

where  A3-  is  the  angle  between  two  adjacent  sample  points  on 
the  same  ring. 

A3-   =  & 

pi 

This  suggests  a  division  of  each  ring  into  an  even  number  of 
patches  NPI  which  will  be 

TTp  . 

NPI   =   2  x  int(-jj-)  (4.3) 

The  function  INT  means  the  integer  value  of  its  argument  so 
that  the  increment  angle  for  the  i   ring  becomes 

A8i  -    m  <4-4) 

Thus,  following  this  procedure,  any  arbitrarily  shaped  body 
of  revolution  can  be  subdivided  into  patches  of  approximately 
equal  surface  area.   Figure  4.2  shows  the  geometry  of  the 
space  sampling  on  the  whole  body. 
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Figure  4.1.   Geometry  of  the  contour  sampling 
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space  sample  point 


Figure    4.2.       Space    sampling   of    the   whole   body 
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B.   NUMERICAL  SOLUTION  OF  THE  INTEGRAL  EQUATIONS 

The  equations  to  be  solved  are  given  by  Eq.  (3.8)  and 
(3.9)  in  previous  sections. 

For  the  numerical  solution  the  currents  will  be  represented 
by  their  space-time  samples: 

J(p,z,t)   =   J(i,n) 

where  i  is  the  i    space  sample  on  the  contour  function  and 
n  is  the  n    time  sample 

i   =   1 ,  .  .  .  ,  NS 
n   =   1 , . . . , NT 

t   =   n-AT 

The  time  being  sampled  at  a  rate  AT.   Similarly  the  incident 
field  is  sampled  and  will  be  denoted  by  its  n-th  time  sample 
at  the  i-th  segment: 

H1(z,t)   =   H^i,!*) 

The  general  procedure  for  the  numerical  procedure  was  already 

mentioned  in  Chapter  II,  and  will  be  reviewed  in  more  detail 

here.   The  incident  field,  H,  is  assumed  to  be  zero  at  all 

times  less  than  some  reference  value  (which  will  be  set  equal 

to  1  At)  .  Hence  the  total  current  on  the  surface  of  the  scatterer 

is  also  zero  prior  to  time  t.  . 

At  the  next  time  interval  tn  =  t  +  At  the  incident  field 

1    o 

reaches  the  scatterer  at  z  =  z,   (first  space  sample)  ,  thus  at 
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time  t,  the  current  will  be  simply  twice  the  tangential  com- 
ponent of  the  incident  field.   As  time  marches  on,  the  incident 
field  covers  more  of  the  body  and  the  current  will  be  due  to 
the  incident  field  plus  the  fields  due  to  currents  at  other 
points  on  the  scatterer  at  earlier  times.   So,  by  marching 
on  in  time,  the  currents  at  each  space  sample,  for  each  time 
increment ,  can  be  computed . 

The  currents  appearing  inside  the  integrals  in  equations 
(3.8)  and  (3.9)  are  the  retarded  currents,  and  are  given  by 

J'   (p,fz\T)   =  J'  ,  (p1  ,z',t-R/C) 
s  ,  cp  s  ,  cp 

The  point  (p',z')  is  represented  by  the  k-th  space  sample 
while  the  time  x  will  be: 

x   =   t  -  R/C   =   n-At  -  R/C   =   gAT. 

In  the  general  case  x  will  not  be  an  integer  multiple  of 
AT,  so  that  the  current  cannot  be  represented  by  any  of  the 
previously  computed  samples.   In  order  to  find  the  value  of 
the  current  at  time  x ,  there  is  a  need  to  interpolate  between 
the  sampled  known,  discrete  values. 

For  a  better  accuracy  use  was  made  of  the  four  point 
Lagrangian  interpolation  formula  using  two  samples  at  each 
side  of  the  value  gAT  if  t,  ,  t_ ,  t  ,  t.  are  the  four  samples 
used  for  interpolation,  such  that: 


t4   =   t   +  At   =   t2  +  2At   =  t1   +    3At 


t4    =    g4  At 
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t3  =  g3  At 
t2  =  g2  At 
t±     =     g1  At 

g-i  /  <3?r    <3r>'    9 a    are  integer  values,  and 
x   =   g  At     (g  real  number) 
such  that 

tL  <  t2  <  T  <  t3  <  t4 

The  interpolation  formula  using  those  four  points  will  be 
given  by:  [17] 

Js,*(T)   "    V-V^-^.*^'1  <4'5) 

n=l  1=1   n   x 

i^n 

This  interpolation  formula  must  be  adjusted  in  two  cases: 

1.  Whenever  t^  or  t.  is  greater  than  the  actual  time 
of  computation  t  (in  order  not  to  interpolate  in 
the  future) . 

2.  Whenever  t.  or  t„  are  less  than  the  beginning  time 
of  computation  t   (due  to  causality) . 

The  time  derivative  of  the  currents  will  be  simply  expressed 
as  the  derivative  with  respect  to  x  of  the  interpolated 
value,  which  will  give  the  following  formula: 
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3 J       ,(t)  4       0         4        T-t 


n=l  i=l      n      l 


3x 

n=± 


Defining: 


Ti    =    g  -  gx 
T2    =    ^  -  g2 

T3    =    g  "  g3 
t4    =    g  -  g4 

The   interpolation   formulas    for   the   current   and   its    time 
derivative   can   be  written   in   their   numerical   form  as    follows 

TxTxT  TxTxT 

,.       .  1         2         3    _        ,.         s  2         3         4 

J       ,    i,g        =      7 J       ,(i  ,g.)     -    -p. 

s,4>  6  S,cp     '^4  6 

T      x  t  .    x  t, 

•J        (i,g,)    +  t J      .  (i,g~) 

s ,  <p         31  2  s  ,  <p         32 

T      x  T      x  T 

"  ^-T1 -Js,<0(i^3)  {4'7) 

and 

3Js,cp(i'g>              T1T2+T2T3   +T1T3    T         ,.  . 

DJs,*      =      3^ =      6AT Js,(P(l'g4) 

T2T3+T3T4+T2T4 
6AT Js,<j)(l'gl) 

T3T1+T4T1+T3T4 

+   2AT Js,^,(l'g2) 

T4T1+T1T2+T2T4 


2AT^ —  Jb,^1'^  (4-8) 
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Now  we  can  write  the  numerical  expression  for  the  currents 
in  equations  (3.8)  and  (3.9)  by  representing  the  continuous 
integration  by  numerical  summation,  and  expressing  the  cur- 
rents in  terms  of  their  space-time  samples. 

For  simplicity,  we  will  write  first  the  expression  for  the 
integral  part  of  equation  (3.8)  and  (3.9)  which  will  be 

denoted  by  J   ,  and  J  , . 
1       cs'       Ctj) 

NP 
NS  I 


Jcs(k'n)   =  W      Z   ^P-^  I 


i=l  m=l        R2 

mj^l  for  I=k  kirn 

V-sI2(m)J  (i,g)+F1CO(m)Js(i,g) 

kirn 

V-sI2(m)DJ  (i,g)+F,CO(m)DJ  (i,g) 
+  4 ^-J: ? ]  (4.9) 


NP. 
NS 


*Wk'n)   -  ^      I      A3TPTAST  I 


i=l  m=l        R,  . 

,,  ,-  i    i       kim 
m^l  for  x=k 


U  -  b  J  " 

[ 


U.Sl'(m) Jg(ifg)+F2CO(m) J  (ifg) 


R,  . 
kim 


U-SI2(m)DJ  (i,g)+F,CO(m)DJ. (i,g) 
+  ^ - 1 2 ]        (4.10) 


where : 


J    (i,g)   are  given  by  equation  (4.7)  and  DJ    (i,g) 
s,(p  by  equation  (4.8)  .  S'^ 

Rkim  =  (4+pl-2PkP.CO(ra)+(z.-zk)2)1/2 
CO(m)   =   Cos [ (m-1) A3- 1 
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Sl(m)   =   sin [ (m-1) Ag. ] 


V  =   Z.  -  z 
1    u 


U   =   V  sin  a,  sin  a .  +  p,  cos  a,  sin  a  .  -p  .  sina  ,  cosa  . 

K       IK       K         11       K       1 

Fn   =   V-CO (m)  -sin  a.  +  (p.  -p . CO (m) ) cos  a. 
1  l     k   l  l 


F_   =   cos  a,  (p,  cos  (m)  -  p  .  )  +  V  sin  a,  CO  (m) 


NS   represents  the  total  number  of  segments 


NP .   represent  the  total  number  of  patches  at  the 
ith  ring. 


The  total  currents  can  now  be  expressed  in  terms  of  J    and 
J  , .  which  yields: 

J_(k,n)   =  -2H1(zv   )  +  J   (k,n)  (4.11) 

s  K,n     cs 

J,(k,n)   =   -2H1(z1  ,n)  *  sin  (a.  )  +J  .  (k,n)       (4.12) 
<p      '  k  k    ccf> 

J   and  J   are  the  values  of  the  current  components  stored  and 
used  for  computation  of  the  backscattered  far  field,  which 
will  be  expressed  in  numerical  form  as  follows: 


s  1    NS 

H  (n)   =   "  4rC   E  [DJ  (i,g)sin  a±+DJ  (i,gJ!ASiPi      (4.13) 
o   i=l    '  "    ? 

where 

z . 

l 
g   =   n  - 


cAt 

and  DJ    (i,g)  are  given  by  equation  (4.8) 

s  /  cp 
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Equations  (4.7)  and  (4.8)  used  for  interpolation  of  the 
currents  and  their  time  derivative  are  given  for  the  general 
case  where: 

1.  g..  >  1  (after  the  incident  field  had  reached  the 
specific  point  i  on  the  scatterer) 

2.  g.  <  n  (so  that  no  interpolation  in  the  future  occurs) . 
In  cases  where  these  two  conditions  are  not  met  the  interpo- 
lation must  be  adjusted  to  include  only  those  points  which 
satisfy  causality  with  no  interpolation  in  the  future. 

The  numerical  formulations  derived  in  this  section  were 
programmed  in  FORTRAN  to  compute  the  currents  and  the  back- 
scattered  field,  due  to  the  excitation  of  a  body  of  revolution 
by  an  impulse  type  field  (or  a  ramp) .   The  input  data  to  the 
program  is  the  function  p(z)  along  with  the  length  of  the 
spatial  segment  AS  and  the  time  increment  AT.   The  program 
is  given  in  Appendix  E. 

C.   NUMERICAL  EXAMPLES 

The  first  example  to  be  considered  is  a  sphere  of  radius 
1  meter.   This  special  example  was  chosen  since  its  result 
appears  in  various  publications  and  comparison  can  hence  be 
made. 

For  the  present  example,  use  is  made  of  the  guidelines  set 
in  Section  (EI. B)  for  the  choice  of  the  incident  field,  the 
sampling  rate  and  the  space  sampling .   The  incident  field  for 
the  case  of  the  impulse  response  is  represented  by  the  gaussian 
shape 
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-A2(t-t    -z/c)2 
max 


2a 


<  t  <  0 


Hx(z,t)   = 


,    otherwise 


where 


A   =   6-108  sec  -1 


which  yields 


AT   =   0.2/C  sec 


The  pulse  is  described  in  figure  4.3, 


Figure  4.3.   Gaussian  shaped  pulse 
incident  on  a  sphere 

The  pulse  reaches  its  maximum  value  at  the  center  of  the  sphere 

t      =   a/C  . 
max       / 
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The  spatial  segment,  with  15  segments 
As   =   ff  =   0-209  m  <  C-AT 

The  number  of  time  increments  observed  was  such  that  it 
covers  6  times  the  transit  time  across  the  whole  sphere  (60-AT). 
The  results  obtained  for  the  backscattered  field  are  given  in 
figure  4.4,  on  which  the  results  of  Bennett  [7]  are  superimposed 
for  comparison. 

Observation  of  the  impulse  response  of  the  sphere  shows 
the  following  results: 

1.  The  first  peak  appears  after  a  transit  time  across 
one  radius  of  the  sphere  (which  is  the  time  where  the  inci- 
dent pulse  reaches  its  peak  value) ,  and  is  due  to  the  return 
from  the  leading  edge  of  the  sphere . 

2.  The  second  peak  appears  at  about  5.2  a/C  and  is  due  to 
the  return  from  the  back  side  of  the  sphere,  after  a  transit 
time  of : 

a  +  I*  +  a   _   5  2- 

c   c   c    ^-z   c 

which  corresponds  to  twice  the  free  space  transit  time  up  to 
the  shadow   boundary  plus  the  transit  time  around  the  back 
of  the  sphere.   This  second  peak  is  the  creeping  wave  return. 

3.  The  response  as  expected,  starts  decaying  after  some 

time  equivalent  to  5-6  times  the  transit  time  across  the  sphere. 

The  ramp  response  for  the  same  sphere  is  given  in  figure 
4.5.   The  first  positive  incursion  of  the  response  waveform 
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radius  sphere 
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extends  to  about  5  times  the  transit  time  across  the  radius 
which  corresponds  to  1.25  round  trips  across  the  target.   The 
impulse  response  of  a  30  cm  diameter  sphere  was  measured  in 
the  NPS  scattering  range  lab.   The  waveform  obtained  is  shown 
in  figure  4.6.   The  response  appears  here  with  a  reversed 
sign  since  vertical  polarization  is  used  (instead  of  the 
horizontal  polarization  considered  for  the  computations) . 
The  waveform  shows  very  good  agreement  with  the  form  obtained 
for  the  2  meter  diameter  sphere. 

The  second  example  considered  was  the  cone-sphere  whose 
geometry  is  shown  in  figure  4.7.   The  target  contour  was 
divided  into  12  segments  while  all  other  parameters  remained 
the  same  as  for  the  sphere.   The  impulse  response  of  the  cone- 
sphere  target  is  given  in  figure  4.8.   The  waveforms  show  two 
positive  peaks,  corresponding  to  the  leading  edge  and  creeping 
wave  returns.   The  separation  between  them  is  equivalent  to  the 
two  way  free  space  transit  time  up  to  the  shadow  boundary 
plus  the  transit  time  across  the  rear  of  the  sphere. 
[(2  xl.3  +i  x0.7)   =   4.7  light-meters]   The  ramp  response 
of  the  same  cone-sphere  target  is  given  in  figure  4.9,  and 
the  remarks  made  on  the  sphere  apply  here  too. 

The  third  example  was  the  capped  cylinder  target  shown  in 
figure  4.10.   The  contour  of  the  target  was  divided  into  12 

segments,  the  sampling  time  was  AT  =  0.025  m/sec  and 

9    -1 
A  =  4 «10   sec   .   The  smoothed  impulse  (gaussian)  response 

of  this  target  is  given  in  figure  4.11.   The  waveform  shows 
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clearly  the  responses  from  the  leading  edge  and  the  rear  of 
the  cylinder,  with  the  same  amplitudes,  while  from  the  flat 
horizontal  region  of  the  cylinder  the  response  is  zero  as 
expected  from  the  analytical  expressions  of  the  field. 

The  ramp  response  for  the  capped  cylinder  is  shown  in 
figure  4.12.   The  first  positive  incursion  of  the  response 
clearly  exhibits  a  shape  comparable  to  the  shape  of  the  target. 
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Figure    4.7. 


Geometry  of    the   cone-sphere 
target    (half-target) 
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Figure    4.8.      Smoothed   impulse    response   of 
the    cone-sphere 
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Figure    4.9.         Ramp   response    of    the    cone-sphere 
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Figure  4.10.   Geometry  of  the  capped  cylinder 
target  (half-target) 
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Figure    4.11. 


Smoothed   impulse   response   of   a 
capped   cylinder    (length    ~    0.33   m) 
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V.   INVERSE  SCATTERING 

A.  PROBLEM  DEFINITION 

Given  the  incident  field  and  the  scattered  field  from  an 
unknown  target,  it  is  desired  in  the  most  general  case  to 
retrieve  the  shape,  size,  orientation,  and,  perhaps,  composi- 
tion of  the  target.   This  is  known  as  inverse  scattering,  or 
target  imaging. 

For  the  present  case,  the  problem  is  restricted  to  a 
metallic  body  of  revolution  illuminated  by  a  field  propagating 
along  the  axis  of  symmetry,  and  the  scattered  field  is  measured 
in  the  backscattered  direction,  as  shown  in  figure  5.1. 

The  solution  of  this  problem  is  approached  here  directly 
in  the  time -domain. 

B.  FORMULATION  OF  THE  PROBLEM 

In  [6] ,  Kennaugh  and  Moffatt  derive  the  relationship  between 
the  physical-optics  impulse  response  of  a  body  and  its  area 
function  A(z),  projected  in  the  plane  orthogonal  to  the  direc- 
tion of  propagation  of  the  field  (z-axis) . 

The  impulse  response  is  given  by: 


FT(t)   =     C   d'A(2) 


Iv  '      87Tr    ,  2 
o   dz 


(5.1) 
z=ct/2 


where  r   is  the  range  from  the  observation  point  to  the 
scatterer,  as  shown  in  figure  5.1.   Integrating  twice  (5.1) 
with  respect  to  time  gives  the  ramp  response  as  follows: 
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Figure  5.1.   Inverse  scattering  geometry 


69 


FR(t)   =   2i7"^A(z) 


(5.2) 
z=ct/2 


Equations  (5.1)  and  (5.2)  were  proven  to  be  valid  for  bodies 
having  a  monotonic  area  function  with  a  limiting  value  at  the 
shadow  boundary.   This  relationship  between  the  ramp  response 
and  A(z)  is  derived  in  Appendix  F  for  the  body  of  revolution. 
The  total  backscattered  field  that  would  be  measured  will 
consist  of  the  physical  optics  term  plus  an  additional  term 
due  to  the  field  induced  interactions  between  the  currents 
at  different  parts  of  the  body.   The  total  current  induced  on 
the  body  is  given  by  equation  (3.1).   The  first  term  of  the 
R.H.S.  of  the  equation  represents  the  physical-optics  current, 
which  we  will  denote  by  J  ,  and  the  second  term  expressed  by 
the  integral  represents  the  additional  current,  which  we  will 
denote  by  Jp.   The  total  current  will  then  be: 

5total(?'t)   =  V?'t}  +  ^C^'t]  (5-3) 

Expressed  in  terms  of  (5.3)  the  total  backscattered  field  for 
a  ramp  excitation  will  be: 


Hl(t»   ■  nhr/   k'V^c'  x^dA'  (5-4) 


=   H^(t)  +  H*{t) 
p        C 


where : 


5 

H  (t)   is  the  physical  optics  ramp  response, 
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s 
Hr(t)  is  the  additional  response, 

and 

g 

H  (t)  is  the  total  ramp  response. 

Using  the  relationship  between  the  physical-optics  response 
and  the  area  function,  (5.4)  can  be  rewritten  as: 

HR(t)   =  2Tcp2{z)    +   HC(t)  '    z  =  ct/2  (5.5) 

o 

which  leads  to 


p2(z)   =   2rQC(H|(t)  -  Hj(t))  (5.6) 


Equation  (5.6)  is  the  basis  for  the  solution  of  the  inverse 
scattering  problem  for  the  body  of  revolution. 

C.   INVERSE  SCATTERING  SOLUTION  PROCEDURE 

Equation  (5.6)  relates  the  contour  function  of  the  body, 
p(z),  to  its  ramp  response.   Solution  of  that  equation  leads 

to  the  shape  of  the  body.   The  only  term  known  in  equation 

s 
(5.6)  is  the  measured  value  of  the  total  ramp  response  H  (t) 

while  Hr (t)  is  not  known. 

Moreover,  all  the  terms  in  the  equation  are  functions  of 

p(z),  hence  there  results  a  non-linear  equation  that  will  be 

solved  by  an  iterative  procedure.   As  a  first  approximation, 

the  unknown  term  Hr(t)  will  be  neglected  and  so  the  first 

estimate  of  the  contour  function  will  be  given  by: 

P;L(z)   -   [2rQC  H^(t)]1/2  (5.7) 
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Once  a  body  is  defined  its  ramp  response  can  be  computed  since 

we  know  the  incident  field.   The  ramp  response  computed  for 

the  first  estimate  of  the  body  will  be  denoted  by  H   (t)  and 

Rl 
is  expressed  as: 

P1(Z) 

HRn(t)  -  -ITc+   Hc\(t)  (5'8) 

1  o       1 

A  new  estimate  of  the  shape  p-(z)  can  be  found  using  (5.6) 
and  (5.7),  and  is  given  by 


p^z)   =   2rQC[H^(t)  -  Hj  (t)]  +  pj(z)  (5.9) 

or  equivelently  by 

p*(z)   =   2rQC[2H|(t)  -  H^  (t)  ]  (5.10) 

The  new  estimate  P2(z)  serves  to  compute  a  new  ramp  response, 

H   (z) ,  from  which  an  updated  estimate  of  p(z)  is  extracted. 

R2 
Following  this  iterative  procedure,  the  k-th  estimate  of  the 

contour  function  will  be  given  by: 


pk(z)   =   pk-l(z)  +  2roC[HR(t)  "  HR    (t)]    t5-11) 

k-1 


or  by: 


k-1 
p£(z)   -   2rnC[k-Hp(t)  -  I    H*     (t)]  (5.12) 

£  =  0   £ 


The  iterations  are  continued  up  to  the  point  where  the  differ- 
ence between  the  true  measured  response  and  the  last  computed 
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response  reaches  a  minimum.   Several  runs  of  the  program  for 
various  bodies  has  shown  that,  as  the  iterations  go  on,  the 
mean  squared  error  is  decreasing,  reaches  a  minimum  then  starts 
oscillating.   In  one  case  observed,  the  error  is  diverging 
after  reaching  the  minimum.   Some  representative  results  are 
shown  in  the  next  section . 

The  normalized  mean  squared  error  will  be  defined  as 
follows: 

(H|(t)  -  H^    (t))2 
£k   =   i ¥ (5-13) 

K  (H|(t)r 

2 
e.  is  the  value  that  must  be  minimized  in  order  to  obtain  the 
k 

most  accurate  representation  of  the  shape  of  the  body.   As 
it  is  defined,  the  error  is  equivalent  to  the  difference  be- 
tween the  successive  shapes,  as  shown  by  equation  (5.11).   Once 

the  difference  H_.(t)  -  H,,    (t)  is  minimized,  the  difference 

K        k-1 
2        2 
p,  (z)  -  p,  1  (z)  is  also  minimized. 

D.   NUMERICAL  SOLUTION  AND  EXAMPLES 

For  the  numerical  solution  we  will  use  equation  (5.12), 
which  is  more  suitable  for  the  algorithm  developed. 

We  will  incorporate  the  following  terminology: 

H_„(nAT)  is  the  time  sample  of  the  measured  ramp  response 

while  H   (nAT)  is  the  time  sample  of  the  computed  ramp  response 

Rk 
for  the  k-th  estimate  of  the  body,  and  P1<(z)  is  the  k-th 

estimate  of  the  contour  function.   In  discrete  numerical 

form,  equation  (5.12)  becomes: 


73 


j  k-1 

Pv(z„)   =   2r  C[k-H_M(nAT)  -  I      H_  (nAT)  ]         (5.14) 
K   n        o      km        £=0   R£ 


where 


z    =   n  C  AT/2 

n  ' 

and  for  I    =  0,  H   (nAT)  =  0.   The  mean  squared  error  is: 

Ro 

N  2 

HHRM(nAT)  -  HR    (nAT)] 
,-2      n=l  k-1  ,_  ._. 

£k   =   N ~ (5'15) 

I      H^(nAT) 
n=l 

where  N  is  the  total  number  of  time  samples  observed  and 
computed . 

The  ramp  response  of  the  body  defined  by  p,  ,  (z  )  will 
be  computed  using  the  integral  equation  solution  to  the  for- 
ward scattering  problem,  as  already  described  in  Chapter  IV. 

The  Fortran  program,  given  in  Appendix  E,  has  been  extended 
with  a  subroutine  named  "SHAPE",  in  order  to  implement  the 
inverse  scattering  algorithm.   In  this  subroutine  the  shape 
of  the  body  is  computed  iteratively,  using  (5.14). 

The  space  samples  given  by  (5.14)  are  defined  at  equal 
increments  on  the  z-axis,  while  the  program  developed  for  the 
computation  of  the  ramp  response  has  been  proven  to  work  best 
when  the  contour  was  divided  into  equal  length  segments .   Con- 
sequently, one  of  the  tasks  of  subroutine  "shape"  is  to  find 
the  space-samples  at  the  points  delimiting  segments  of 


74 


approximately  equal  length.   This  is  done  using  a  "cut  and 
try"  method  which  converges  to  the  desired  length  of  segment. 

The  iterative  approach  for  solving  the  inversion  equation 
(5.14)  was  tested  in  a  simulative  way  for  the  three  types  of 
bodies  already  considered  for  the  forward  scattering  solution 
in  the  previous  chapter. 

The  ramp  response  of  the  body  is  first  computed  and  is 
used  to  simulate  the  response  that  would  have  been  measured 
(without  noise  or  error) .   Then  the  shape  of  the  body  is  re- 
trieved from  that  ramp  response. 

The  shape  of  the  body  is  retrieved  from  the  first  positive 
incursion  of  the  ramp  response,  which  as  it  will  be  seen 
extends  beyond  the  length  of  the  target  for  the  shapes 
considered. 

The  first  example  was  the  1  meter  radius  sphere .   Figure 
5.2  shows  the  1-st,  3-rd,  6-th  estimates  of  the  contour  func- 
tion of  the  sphere.   The  6-th  estimate  gives  the  least  error 
and  is  shown  superimposed  on  the  original  shape  in  figure  5.3. 

The  second  example  was  the  cone-sphere,  used  in  two  con- 
figurations.  The  first  one  when  the  contour  of  the  body  was 
originally  cut  to  12  segments.   Figure  5.4  shows  the  successive 
estimates  (1,3,6)  of  the  contour,  and  figure  5.5  shows  the 
minimum  error  estimate  (estimate  6) .   The  second  version, 
when  the  body  was  cut  to  11  segments ,  gave  a  much  more  accurate 
result,  as  shown  in  figures  5.6  and  5.7. 
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The  last  example  was  the  capped-cylinder .   Figure  5.8 
shows  the  successive  estimates  of  the  shape  (estimates  1,3,4), 
and  figure  5.9  shows  the  4-th  estimate  which  gives  the  minimum 
error. 

For  the  case  of  the  capped  cylinder,  after  reaching  the 
minimum  error  at  the  4-th  iteration,  the  imaged  shaped  of 
the  body  started  to  diverge  considerably,  and  appeared  more 
shortened  than  the  actual  shape. 

Common  to  all  examples,  the  error  appears  mostly  at  the 
shadowed  region  of  the  body,  but  in  general  the  shape  is 
retrieved  with  quite  a  good  fidelity,  which  proves  that  the 
iterative  direct  time  domain  solution  is  working,  subject  to 
certain  caveats  which  will  be  discussed  in  the  conclusions. 

E.   INFLUENCE  OF  NOISE 

In  the  examples  considered  in  the  previous  section  the 
only  noise  appearing  in  the  ramp  response  is  the  numerical 
noise  generated  by  sampling  and  computational  round-off.   In 
this  section  the  measurement  noise  will  be  tested  with  a  ramp 
response  to  which  noise  is  added.   In  this  simulation,  each 
noise  sample  is  obtained  by  the  sum  of  a  large  number  of  zero- 
mean,  uniformly  distributed  random  numbers,  which,  by  use  of 
the  central  limit  theorem,  approximate  the  statistics  of 
gaussian  noise.   The  noise  samples  are  normalized  in  order  to 
remain  within  a  certain  voltage  signal  to  noise  ratio,  and 
then  are  added  to  the  computed  ramp  response. 

The  body  considered  for  the  noise  test  was  the  sphere, 
and  the  signal  to  noise  ratios  tested  were: 
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V 
2  0  log  ^^ —  =   10,  2  0,  30   and   4  0  dB  . 
noise 

The  sphere  ramp  responses  with  the  noise  added  are  shown  in 
figures  5.10,  5.11,  5.12  and  5.13  for  the  successive  signal 
to  noise  ratios. 

For  the  10  dB  S/N  case,  the  appearance  of  spikes  and  dis- 
continuities on  the  ramp  response  causes  the  shape  of  the 
body  to  diverge  considerably  since  the  first  iteration,  and 
the  shape  cannot  be  recovered . 

For  the  2  0  dB  s/N  case,  the  computed  shape  starts  to  con- 
verge, with  the  minimum  error  shape  obtained  at  the  third 
iteration.   Then  the  imaged  shape  proceeds  to  diverge  and  to 
shorten. 

Figure  5.14  shows  the  successive  estimates  of  the  shape 
up  to  the  3-rd,  and  figure  5.15  shows  the  sixth  estimate,  which 
is  clearly  different  of  the  actual  shape. 

For  the  30  dB  S/N  case,  the  minimum  error  shape  was  ob- 
tained at  the  5-th  iteration,  as  shown  by  figure  5.16.   Then 
again,  the  computed  shape  diverges  as  shown  in  figure  5.17. 
Common  to  all  cases  of  divergence,  the  wrong  estimates  start 
as  soon  as  the  last  computed  shape  exhibits  a  strong  discon- 
tinuity in  its  contour.   A  very  important  conclusion  from  the 
noise  analysis  is  that  before  starting  the  inverse  scattering 
procedure  the  ramp  response  used  must  be  smoothed  so  that  no 
spikes  or  sudden  discontinuities  appear  on  it. 

We  should  normally  expect  that  the  true  measured  response 
will  have  a  20  dB  s/N  ratio,  so  that  the  problem  of  divergence 
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Figure  5.10.   Ramp  response  of  the  sphere  with 
noise  added  ( S/N  =  10  dB) 
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Figure  5.11. 
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Figure  5.12.   Ramp  response  of  the  sphere  with 
noise  added  (S/N  =  30  dB) 
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Figure  5.13.   Ramp  response  of  the  sphere  with 
noise  added  (S/N  =40  dB) 
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must  be  investigated.   In  any  case,  before  diverging  the 
shape  that  gives  the  minimum  error  is  quite  a  good  approxi- 
mation of  the  true  shape. 

An  additional  test  with  40  dB  S/N  was  made,  and  no  signi- 
ficant deviation,  from  the  case  without  noise  added,  was 
observed.   Figure  5.18  shows  the  first  and  6-th  estimates  of 
the  shape.   The  minimum  error  shape  was  the  6-th  estimate. 
Up  to  the  12-th  iteration  the  shape  of  the  body  remained 
stable,  with  minor  variations. 
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VI.   SUMMARY 

A.   SUMMARY  OF  RESULTS 

The  inverse  scattering  method  for  radar  target  imaging 
was  the  main  topic  addressed  in  this  work.   The  inverse 
scattering  algorithm  developed  uses  the  time  domain  solution 
of  the  integral  equation.   The  class  of  targets  considered 
included  the  conducting  bodies  of  revolution  whose  important 
symmetry  property  reduces  considerably  the  computation  time 
and  the  computer  storage  needed  for  the  solution  of  the  in- 
verse scattering  problem. 

The  use  of  the  ramp  response  for  retrieving  the  shape  of 
the  target  proved  to  be  a  very  efficient  tool,  and  the  exam- 
ples considered  gave  generally  good  results,  with  initial 
convergence  of  the  imaged  shape.   For  the  case  of  the  long 
cylinder  and  the  sphere  with  low  signal  to  noise  ratios  there 
appeared  an  eventual  divergence  of  the  imaged  shape  after 
many  iterations.   This  is  thought  to  be  due  to  accumulated 
numerical  roundoff  errors  in  repeatedly  solving  the  integral 
equations.   In  general  it  has  been  shown  that  the  use  of 
the  ramp  response  by  itself,  without  any  further  computation, 
gives  a  close  approximation  to  the  shape  of  the  body  but  up 
to  its  shadow  boundary  which  corresponds  to  the  maximum  of  the 
area  function  A(z)  along  the  line  of  sight.   Beyond  that  boundary 
the  body  is,  at  the  first  observation,  elongated  and  the 
iterative  algorithm  is  needed  to  compensate  for  the  errors 
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induced  by  the  wake  or  "creeping  wave"  which  alters  the 
target  response. 

As  the  iterative  computation  process  goes  on,  some  numeri- 
cal instability  that  is  mainly  due  to  a  discontinuity  in  the 
computed  shape  might  arise.   However,  before  diverging,  the 
shape  of  the  body  always  reaches  an  estimate  which  has  a 
minimum  error  relative  to  the  true  shape,  and  is  a  good 
approximation  of  that  shape,  so  that  the  minimum  error  cri- 
teria set  for  stopping  the  iterative  process  appear  to  be  an 
attractive  way  to  overcome  the  instability  phenomena.   As  a 
continuation  of  this  work  it  is  suggested  to  carry  out  a 
detailed  analysis  of  the  error  and  numerical  instabilities 
that  have  been  observed  while  testing  the  iterative  algorithm. 

B .   RECOMMENDATIONS 

The  algorithm  developed  for  target  imaging  by  transient 
inverse  scattering,  has  been  tested  numerically  by  simulating 
additive  physical  noise.   It  must  now  be  evaluated  in  a  real 
laboratory  environment  where  the  true  measured  response  is 
used.   The  transient  scattering  range  at  NPS  should  be  used 
to  test  and  finalize  the  algorithm  for  the  class  of  targets 
addressed  in  this  thesis. 

The  first  step  to  be  done  is  to  obtain  the  measured  ramp 
response  of  the  target.   The  range  uses  an  impulse  generator 
which  exhibits  a  gaussian  pulse  shape  characteristic.   The 
true  ramp  response  can  be  obtained  by  using  either  a  time 
domain  convolution  or  a  frequency  domain  computation  and  I.F.T. 
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The  true  ramp  response  can  be  obtained,  in  principle,  from 
the  gaussian  response  of  the  target  using  the  following 
relationship: 

H   U) 

CO 


where : 


HR(oo)   is  the  Fourier  transform  of  the  true  ramp 
response  of  the  target, 

H   (w)   is  the  gaussian  impulse  response  of  the 
target, 

G(co)   is  the  Fourier  transform  of  the  transmitted 
gaussian  impulse, 


and  the  ramp  response  will  be: 


HR(t)   =   F_1[HR((JJ)]  . 


Very  special  care  must  be  given  to  the  choice  of  space 
sampling  and  time  sampling  intervals.   The  time  sampling  is 
dictated  by  the  sampling  rate  of  the  oscilloscope,  and  the 
space  sampling  must  be  done  at  a  high  enough  rate  such  that  in 
every  case  the  distance  R  between  the  closest  spatial  samples 
is  related  to  the  time  sampling   T  by: 

R/C   >   AT  . 

The  next  step  is  to  translate  the  FORTRAN  program  given  in 
Appendix  E  to  fit  the  computer  of  the  lab.   There  is  a  need 
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to  reduce  to  a  minimum  the  computation  time  since  the 
laboratory  microcomputer  is  very  slow  compared  to  a  large 
mainframe.   The  use  of  some  additional  symmetries  of  the 
body  can  certainly  help  in  this  task. 

Finally,  it  is  necessary  to  minimize  the  noise  in  measure- 
ments, and  smooth  the  response  to  be  used  for  shape  computation, 

C.   AREAS  FOR  FURTHER  STUDY 

The  ultimate  objective  of  this  research  program  is  to 
aid  in  determining  the  feasibility  of  developing  future  radar 
systems  having  the  innate  ability  to  classify  targets  using 
their  transient  scattering  response. 

In  the  most  general  case  the  target  will  be  of  arbitrary 
shape  and  the  use  of  the  iterative  method  described  in  this 
thesis  would  be  highly  impractical  and  prohibitive  for  com- 
plex targets,  due  to  two  principal  factors.   The  first  is  the 
time  needed  to  compute  the  response  of  a  complex  body,  and 
the  second  is  that  the  area  function  obtained  from  the  ramp 
response  is  not  sufficient  to  determine  in  what  way  the 
response  is  to  be  computed. 

The  best  technique  to  deal  with  such  targets  would  involve 
a  way  that  does  not  require  the  computation  of  the  integral 
equation,  and  instead  will  use  only  the  measured  impulse 
response  (or  ramp  response)  to  discriminate  between  targets. 
The  first  solution  involves  measuring  the  ramp  response  at 
multiple  look  angles.   The  response  waveform  will  be  used  up 
to  the  shadow  boundary  to  find  the  area  function  along  each 
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line  of  sight.   The  target  image  can  be  approximately  recon- 
structed from  the  projections  of  the  area  functions.   Das 
and  Boerner  [18],  suggested  a  very  interesting  approach  for 
target  shape  estimation  using  these  projections.   This  method 
seems  attractive  and  is  felt  that  it  could  lead  to  a  good 
solution  for  target  imaging.   Its  main  disadvantage  is  the 
need  for  multiple  look  angles,  requiring  either  interfaced 
multiple  radars,  or  long-time  dwelling  with  one  radar. 

An  alternative  approach  is  to  classify  targets  via  the 
complex  poles  and  residues  of  their  time-domain  signatures. 
Once  the  poles  and  residues  are  extracted  they  could  be  com- 
pared with  stored  library  data  corresponding  to  known  targets, 
and  thus  would  provide  a  classification  of  the  target  under 
test.   This  catalogue  approach  is  also  being  pursued  using 
the  NPS  transient  range,  as  well  as  through  other  independent 
efforts,  and  may  prove  to  be  the  most  practically  implementable 
technique  available. 

Another  area  which  could  form  the  subject  of  future  re- 
search includes  the  time  correlation  of  the  shape  of  the 
target  with  its  motion. 
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APFENDIX  A 
DERIVATION  OF  THE  SCALAR  INTEGRAL  EQUATIONS 

In  Chapter  III,  equations  (3.  3)  and  (3 . 4)  involve  some  vector 
multiplications.   The  first  appearing  in  the  integrals  is: 


a  xj'  xr  =   (a  xa1  xr).j'  +  (a  x  a !  xr)j'         (a1) 
n  n   s      s     n   <|>     $ 


We  note  the  following  identities: 


a  xa'  xr  =   a'  (a  -R)  -  R(a  -a1)  (A2) 

n   s         s   n         n   s 


/\  s* 


a  xa'  xr  =   a!(a  .R)  -  R(a  -a!)  (A3) 

n    (£         <p      n         n   cp 


/\  /\      /■* 


Expressing  R,  a1  and  a!  in  terms  of  the  (a  ,a,,a  )  coordinates 
c  ^      s       <p  s   tj)   n 

system  gives. 


R      =      aA (R-a, )    +   a    (R-a    ) 
tp  <p  s  s 

a!      =      a. (a!- a.)    +   a    (a'-a)  (A4) 

<p  <p      <p      <J>  s      <p      s 

a'      =      a. (a'-a    )    +   a    (a'«a    ) 
s  cf>      s      s  s      s      s 

For  the  body  of  revolutions  we  note  the  following  relationships 

a'-a,   =   cos  3 

a!-a    =   -  sin  a -sin  3 
(j)   s 

a!»a    =   -  cos  a-sin  3 
(j)   n 

a'*a    =   sin  a  sin  a 'cos  3  +  cos  a  cos  a' 
s   s 
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a1 -a ,   =   sin  a '  sin  3 
s   <p 

a1 -a   =   cos  a  sin  a 'cos  3  -  sin  a  cos  a' 
s   n 


where 


3  =  <j>'  -  <p 


/\  S\  /N 


The  vectors  R,  a  ,  a  ,  aA  are  given  in  cartesian  coordinates 

n        s        <j> 

by: 


R     =      a    (p  cos  4>  -  p  '  cos  cj) ' )    +   a    (p  sin  4)  -  p  '  sin  <j> ' ) 
x  y 

-  a    (z*  -z) 

z 

a        =      a    (cos  a  -cos    d>)    +   a    (cos  a  «sin   d>)    -   a      sin   a 
n  x  Y  z 

a        =      a    (sin  a  -cos    d>)    +   a    (sin  a  »sin  d> )    +  a      cos   a 
s  x  T  y  z 

a,      =      -a    (sin   a)    +   a    -cos   d> 
<j>  x  y 

Carrying  out  the  dot  products  gives: 


R«a   =  cos  a(p  -p'cos  3)  +  (z1  -z)sin  a 

R»a    =   sin  a(p  -p'cos  3)  -  (z'  -z)cos  a  (A5) 

R-a,   =   -  p'  sin  8 

Using  the  relationships  given  by  the  sets  of  equations  (A4) 
and  (A5) ,  and  simplifying,  equations  (A2)  and  (A3)  become: 
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/\  /\ 


a     xa'   xR      =      a, {sin    3  [  (z  ' -z)  sinasina  '  +pcosasina' 
n        s  <p 

-   p'cosa'sina] }  +a    { ( p-p ' cos3) cosa ' 

+    (z '-z) cosgsina ' }  (A6) 


a     x  a '  x r  =   a  ,  [  (z '-z) sinacos  3+  cosa (pcos3-P ' ) ] 


+  a  [-sing(z'-z)  ]  (A7) 


Defining : 

V      =  z1    -    z 

F.      =  (z  '  - z) cos  3sin   a'    +   cos   a '  (p-p 'cos    3) 

U     =  (z  '  -  z)  sin  a  sin   a  *    +  p  cos  a  sin  a  '  -  p  '  cos  a  '  sin   a 


F?      =      cos   a(pcos  3 -p  ' )    +   sin   a   cos    3(z'  -  z) 


produces    from    (A6)    and    (A7) 


a     xa'   xr     =      a,    U  sin  3    +   a    -F,  (A8) 

n         S  q>  si 


/N  s> 


a     x a '   x R      =      a,    F_    -    a      V-sin    3  (A9) 

n        (p  <p      2  s 


Note    that   V,    F..  ,    F~    are   even    functions   of    3    and   U   is    inde- 
pendent of   3.      Next  we   have   to   carry   the   second   vector  multi- 
plication: 


a     xH        =      -a      H      sin    d>    -    a,    H      cos    <p    sin    a  (A10) 

n  s  y  cp  v 
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where  H   is  given  in  cartesian  coordinates  by: 
H    =   ax-H 

After  substituting  equations  (A8) ,  (A9),  and  (A10)  into 
equations  (3.3)  and  (3.4)  we  obtain  the  scalar  expressions 
for  the  components  of  the  current  in  the  s  and  <p  directions 
as  follows: 

J   =  -2^31114,  +  ±   Jl[|  +  1  l-iEj^-J^Vsin  6]dA-   (All) 

A  R  Y 

J        =      -2H1cos  (J)  sin  a   +  ^     jl[l+l  fr1 

y  A   R 

•[J*    U   sin    6   +  J'F„]dA'  (A12) 

S  (p     Z 
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APPENDIX  B 

DEMONSTRATION  OF  THE  SINUSOIDAL  VARIATION 
OF  THE  CURRENTS  WITH  <j> 

Maxwell's  first  equation: 

■+ 
VxH   =   £  E  +  J  (Bl) 

In  the  case  of  interest,  due  to  symmetry,  the  field  H  is 
independent  of  tj>  and  is  given  by 


H   =   a   H  (p,z,t) 
x   o 

or  in  cylindrical  coordinates 

H   =  H   cos  d>  a   -  H   sin  d>  a, 

o  p      O  <j> 

+13                                   A           ^Ho 
VxH      =       [—  -r — (p  H     )  ]sin  d>    a      +   -* cos  d>  a, 

p  3z  K  o       Y   p    3z      T(j) 

■  ,     3Ho  (B2) 

-  a   sin  d)  •-* — 
z        3p 

The  currents  on  the  body  expressed  in  cylindrical  coordinates 

are: 


J  =   a   sin  a  J  +  a,  J,  +  a   cos  a  J  (B3) 

p        s     (b   <b     z        s 

The  electric  field: 


E      =      a      E        =      a      sin    <b    E      +   ax    cos    <b    E  (B4) 

yo  p  To<j)  Yo 

Substituting  equations   (B2) ,  (B3)  and  (B4)  in  Maxwell's 
equation  (Bl)  and  equating  the  terms  in  the  same  direction  gives. 
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1  9 
sin  orJ    =   [—  ir—  (pH  )  -  E  ]  sin  <j) 
s       p  9z    o      o 

The  term  inside  the  brackets  is  independent  of  <j>  so  that  we 
can  write 


sin  a  J  =  J   (p,z,t)  sin  A  (B5) 

so. 

Similarly 

9H 
J,   =  (-^   -  E  )cos  <t> 

J,   =   J   (p,z,t)  cos  <p  (B6) 

<p  o2 

and  finally 

9H 
o 
J   cos  a   =   -  sxn  a 


9po 


=   J   (p,z,t)  sin  <J)  (B7) 

°3 

Equations  (B5) ,  (B6) ,  (B7)  show  clearly  that  the  currents 

J  and  J,  vary  sinusoidally  with  the  angle  <p    so  that  we  can 

generalize : 


J  (p,z,t,<J>)   =   J  (p,z,t)  sin  <j> 


J   (pfZ,tf(J>)    =   J  (p,Z,t)   COS  <J) 
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APPENDIX  C 
SIMPLIFICATION  OF  THE  SCALAR  INTEGRAL  EQUATIONS 

In  Chapter  III  the  scalar  equations  for  the  currents 
components  were  given  by  equations  (3.4)  and  (3.5).   In  these 
equations  we  had: 

dA'   =   p'  ds'  dB 
and 


Js   =   J   sin  d> 
s 


Jd>   =   J,  cos  <b 


Substituting  in  equation  (3.4)  gives 

Js  sin  4>  =  -2H1sin  <\>   +  j^     /  P'ds'   /   -yf^  87  +  R] 


2tt 
0         R* 


•[J'    sin<J>'F,  -J'    sin  6    Vcos<t>']dB  (CI) 

o  X      CD 


As  noted  in  Appendix  A,  the  expressions  V,  F.  and  F   appearing 
in  the  integral  are  even  and  periodic  functions  of  S.   We 
recall  that 

?     2  2  1/2 

R   =   (pz  +  p'   -  2ppcos  B  +  (z  '  -  z)  )  x/z 

and 

x   =   t  -  R/C 
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so  that  R,  J'(p',z,x)  and  J'(p',z',t)  are  also  even  and 
periodic  functions  of  3-   Finally  U  is  independent  of  g. 

We  will  concentrate  on  the  integral  over  3  in  order  to 
reduce  it  to  its  simplest  expression.   Since 

<J>'   =  ij)  +  B 

we  can  write 

sin   <J>  *  =      sin   $    cos    3   +   sin    6   cos    <j> 

cos    <j> '  =      cos    <j>   cos    3    -   sin    <j>   sin    3 

and  substitute  this  into  the  integral  over  3  in 
equation  (Cl) .  The  integral  over  3  for  a  complete  period,  of 
all  terms  that  contain  sin  3  will  be  zero  since  those  terms 
are  periodic  and  odd  functions  of  3.   All  the  other  terms  will 
produce  a  non-zero  integral  and  will  be  retained. 

After  this  simplification,  equation  (Cl)  will  become: 

1  s 

J      sin    cj>      =      -2H      sin    d>    +    sin    d>    -=—      (    p'ds' 
s  2ir      J 


2tt 
0    R 


'   "Ml:r  +  R'3  [JsFl  COS  6  +VJlsin23]dB    (C2) 


which  after  repressing  the  common  term  sin  cj> ,  leads  to  the 
final  equation: 

Js   ■   -2Hi  +  ^A'ds'  J%hlr  +  b 

V  (J  K 

•  [J'F.  cos  3  +  VJ!sin23]d3.  (C3) 

si  <p 
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Following  the  same  reasoning  for  equation  (3.5),  and  repressing 
the  common  term  cos  cf>  leads  to  the  following  final  equation: 

•[J1  U  sin26  +J!F_cosg]d3  (C4) 

s  (p  2 

In  both  equations  (C4)  and  (C3)  the  currents  are: 


J   ,   =   J   ,  (p ,  z, t) 


J'     =   J»   (p1  ,z,  ) 

Sf<J)        S,(j> 

t  being  the  retardated  time  t  -  R/C,  and  the  incident  field 
is: 


H1   =   H1(z,t) 
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SYSIN  DO  DATA 

APPENDIX  D 

C 

c 
c 
c 

Q  ************************************************** 

C         *   TRANSIENT  RESPONSE  OF  A  STRAIGHT  THIN  WIRE     * 
C         *  * 

q         **********  ***************************  ************* 

C 

C 

c 
c 

C        PURPOSE: 
C        ******** 

C         THIS  PROGRAM  COMPUTES  THE  TRANSIENT  TIME  DOMAIN 

C  RESPONSE  OF  A  STRAIGHT  THIN  WIRE. 

C 

C         THE  PROGRAM  COMPUTES  THE  CURRENTS  INDUCED  ON  THE 

C  WIRE  SY  AN  INCIDENT  IMPULSE  FIELD, AND  THE  SCATTERED 

C  FIELD,  FOR  ANY  ANGLES  OF  INCIDENCE  AND  SCATTERING. 

C 

c 

c 

C       USAGE: 
r  ****** 

C         THE  PROGRAM  CAN  BE  USED  AFTER  PROPER  DEFINITION 

C      OF  THE  INITIAL  PARAMETERS  OF  THE  WIRE, AND  DEFINITION 

C      OF  THE  DESIRED  ANGLES  OF  INCIDENCE  AND  REFRACTION. 

C 

C         THE  WIRE  IS  ALIGNED  WITH  THE  Z-AXIS  STARTING  AT  Z=0 

C      THE  ANGLES  ARE  MEASURED  WITH  RESPECT  TO  THE  Z-AXIS 

C 

C 

C 

C         PARAMETER  DEFINITION 

C         ******************** 

C 

c 

C         8ETA<400,40)  :ARRAY  DEFINING  THE  CURRENTS  FOR  UP 
C  TO  40  SEGMENTSAND  400  TIME  INCREMENTS. 

C         GAMMA!  400, 40)  : ARRAY  DEFINING  THECHARGES  FOR  UP 

C  TO  40  SEGMENTS  AND  400  TIME  INCREMENTS. 

C 

C         NS:  NUMBER  OF  SEGMENTS  IN  WHICH  THE  WIRE  IS  DIVIDED 

C 

C         SL:  TOTAL  LENGTH  OF  THE  WIRE  IN  METERS. 

C 

C         A:  RADIUS  OF  THE  WIRE  IN  METERS. 

C 

C  ITS:    NUM3ER   OF    TIME    INCREMENTS. 

C  ITS    MUST    BE    A    MULTIPLE    OF    80    FOR    PLOTTING 

C  CONVENIENCE 

C 

C         TI:  ANGLE  OF  INCIDENCE  IN  DEGREES. 

C 

C         TS:  ANGLE  OF  REFRACTION  IN  DEGREES. 

C 

C         DZ:    LENGTH  OF  A  SEGMENT  (L/NS) 

C 

C  DT:    TIME    INCREMENT    IN    SECONDS    (DT=DZ/C) 

C 

C  AA:    GAUSSIAN    IMPULSE    PARAMETER. IN    SEC-1. 

C 
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C         TMAX:  TIME  DELAY  OF  THE  IMPULSE  FROM  THE  FIRST  TIME 

C  OF  INCIDENCE  TO  ITS  PEAK. ( TYPICALLY  DT*NS/2) 

C 

C         EI:  INCIDENT  FIELD  (GAUSSIAN  PULSE). 

C 

C  ESKOO):    ARRAY    CONTAINING    THE    COMPUTED    SCATTERED 

C  FIELD    UP    TO    400    TIME    INCREMENTS. 

C 

C  DELAZ    :    TIME    DERIVATIVE      OF    THE    MAGNETIC    POTENTIAL. 

C 

C         DELPHI:  TIME  DERIVATIVE   OF  THE  ELECTRIC  POTENTIAL. 

C 

C  C:    SPEED    OF   LIGHT. 

C 

C  Cl:    CONSTANT    (MU/4*PI). 

C 

C  C2:    CONSTANT    ( 1 /4*PI *EPS ) . 

C 

C  X(30J,Y(30):C00RDINATES    REQUIRED    FOR    SUBROUTINE 

C  'PLOTP1. 

C 

C 

C  REMARKS: 

C         THE  ANGLES  OF  INCIDENCE  AND  REFRACTION  DESIRED  ARE 

C  ENTERED  IN  DATA  CARDS  ACCORDING  TO  FORMAT  2F6.1. 

C  EACH  DATA  CARD  CONTAIN  ONE  ANGLE  OF  INCIDENCE 

C  AND  ONE  ANGLE  OF  REFRACTION. 

C  THE  COMPUTATION  IS  DONE  FOR  AS  MANY  DATA  IMPUTS  AS 

C  THERE  ARE  DATA  CARDS. THE  LAST  DATA  CARD  MUST  CONTAIN 

C  THE  DATA  NUMBER  1.0  IN  ORDER  TO  TERMINATE  THE  PROGRAM. 

C 

C 

C 

C         SUBROUTINES  USED: 

C         THE  PROGRAM  USES  A  LIBRARY  SUBROUTINE  NAMED  PLOTP 

C      FOR  THE  PLOT  OF  CURRENTS  AND  FIELD  AS  A  FUNCTION 

C      OF  TIME. 

C 

C 

c 

C         METHOD  OF  SOLUTION: 

C         THIS  IS  THE  SOLUTION  OF  THE  E.F.I.E.  EQUATION 

C  APPLIED  TO  A  THIN  WIRE , ACCORD ING  TO  THE  METHOD 

C  DESCRIBED  BY  HARRINGTON  AND  SAYRE  IN  'TIME  DOMAIN 

C  RADIATION  AND  SCATTERING  BY  THIN  WIRES  ■  (IN  »APPL. 

C  SCI. RES.  26*  SEPT  1972) 

C 

C 

c 
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c 
c 

Q  ***************** 

C       *  MAIN  PROGRAM   * 

C       *  * 

C       ***************** 

C 

c 
c 

DIMENSION   BETA!  <*00, 40),  GAMMA  (400, 40),  ESC  400) 

DIMENSION  XC400),Y(400) 
C 
C 

C  CONSTANTS    INITIALIZATION 

C 

C=3.0E+8 

Cl=1.0E-7 

C2=9.0E+9 

PI=3. 1^1592654 
C 
C 

C  PARAMETERS    IMI TI ALIZATION 

C 

C  COMPUTATION    PARAMETERS 

C 

NS=40 

ITS=400 

NSl=NS-l 
C 

C  WIRE    PARAMETERS 

C 

SL=1.04 

A=7.75    E-4 

DZ=SL/NS 

0T=0Z/C 
C 
C 

C  GAUSSIAN    IMPULSE    PARAMETERS 

C 

AA=0.67    E+9 

TMAX=OT*L8.0 
C 

10  CONTINUE 

C 

C  REAOING  OF    DESIRED   ANGLE    OF    INCIDENCE    TI 

C  AND    ANGLE    OF    REFRACTION    TS. 

C 

READC5.9510)     TI,TS 
9510       F0RMAT(2F6.  1) 

IF(TI.EQ.l.O)    GO    TO    999 
C 

C         WRITE  HEADING 
C 

WRITE(6,9601)     TI,TS 
9601       FORMATl 1H1,10X, 'CURRENTS    INOUCED    WHEN    FIELD    INCIDENT', 
2*AT    ANGLE* ,1F6.1,//, 10X, 'FAR    FIELD    SCATTERED    TOWARDS*, 
3 'ANGLE* , 1F6. 1, ////, 8X, *TI ME *,15X,* CURRENT    (AMP)*,///) 
C 

C  TRANSLATION    FROM    DEGREES    TO    RADIANS. 

C 

TIR=PI*TI/180.0 

TSR=PI*TS/180.0 

AI=GOS<TIR) 

AR=COS(TSR) 
C 
C 
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C         START  OF  CURRENTS  COMPUTATION 

C 

c 

Sl=0.5 

S2=-0.5 

S3=(A/DZ)**2 

FO=ALOG({  S1  +  SQRT(  Sl**2+  S3 1  )  /(  S2+SQRT(52**2+S3 ) )  ) 
C 

C  INITIAL    CONDITIONS    SETTING 

C 

DO    100   L=1,NS1 

L1  =  L-L 

IF(TI.LE.90.0)    SO    TO    120 

L2=L-NS 

GO  TO  130 
120    CONTINUE 

L2=L 
130    CONTINUE 

EI=  EXP((-AA**2)*(DT*( 1-L2*AI )-TMAX)**2)  *SIN (TIR) *3.0 

BETA(1,L)=EI/(2*C1*F0/DT) 

IF<L.EQ.l)   GO  TO  110 

GAMMA U,L)=-( 1. 0/C )*( BETA ( 1,L) -BETA! 1, LI) I 

GO  TO  100 
C 

C         BOUNDARY  CONDITIONS  SETTING 
C 
110    CONTINUE 

GAMMA<1,1  )=-BETA  (1,1   ) /C 
100    CONTINUE 

GAMMA(1,NS)=  3ETA  (l,NSl)/C 
C 

C         START  ALL  CURRENTS  COMPUTATION 
C         AT   EACH  TIME  INTERVAL, AND  FOR  ALL  SEGMENTS 


C 

c 
c 


DO  200  J=2,ITS 
DO  300  L=1,NS1 


DELPHI      =0.0 

DELAZ       =0.0 

DO  310  K=1,NS1 

K1=IA8S(L-K) 

K2=J-K1 

Sl=Kl+0.5 

S2=Kl-0.5 

S3=<A/DZ)**2 

F  =ALOG((  SH-SQRT(  Sl**2+S3)l/(  S2+SQRTC  S2**2*S3  ) )  ) 

IF(Kl.EQ.O)  GO  TO  320 

IF(K2.LE.O)  GO  TO  310 

DELAZ  =DELAZ  +BETA( K2, K) *F 
320    CONTINUE 

K3=K2-1 

IF(K3-LE.0)  GO  TO  310 

DELAZ  =DELAZ  -BETA( K3, K ) *F 

310         CONTINUE 

DELAZ  =DELAZ  *C1 

DO    330    K=1,NS 

K2=IABS    (L-K) 

K1=IABS    (L-K+l) 

K3=J-1-K1 

K4=J-1-K2 

IF(K3.LE.0I    GO   TO    340 

Sl=Kl+0.5 

S2=Kl-0.5 

F    =  AL0G((SH-SQRT(Sl**2+S3))/(  S2+SQRT(  S2**2+S3  ) )  ) 

DELPHI  =DELPHI  +GAMMA (K3 , K) *F 
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3^0    CONTINUE 

IFCK4.LE.0)  GO  TO  330 

Sl=K2+0.5 

F  =ALOG((  SH-SQRT<  Sl**2+ S3 ) ) /( S2+SQRT( S2**2+S3 ) ) ) 
DELPHI      =DELPHI      -GAMMA <K4, K) *F 

330  CONTINUE 

DELPHI      =OELPHI      *C2 

IFCTI.LE. 90.01  GO  TO   331 
C 

C         FOR  ANGLES  OF  INCIDENCE  MORE  THAN  90  DEGREES 
C 

L2=L-NS 

GO  TO  332 
C 

C         FOR  ANGLES  OF  INCIDENCE  LESS  THAN  90  DEGREES 
C 

331  CONTINUE 
L2  =  L 

332  CONTINUE 

IF  <J.GT.50  )  GO  TO  333 

EI=  EXPU-AA**2)*(  J*DT*(1-L2*AI  )-TMAX)**2>  *SIN(TIR)*3 

GO  TO  33^ 

333  CONTINUE 
EI=0.0 

334  CONTINUE 

BETA<J,L)=( EI-DELAZ  /DT-DELPHI  /DZ ) /( C1*F0/DT ) 

J1=J-1 

L1=L-1 

IFtL.GT.l J    GO    TO    350 

GAMMA(J,L)=    GAMMA(J1,L     )-BETA(J,L       )/C 

GO    TO    300 
350         CONTINUE 

GAMMA(J»L)=    GAMMA(J1,L    )-BETA(J»L       l/C    +BET A( J, LI )/C 
300         CONTINUE 

GAMMA( J,NS)=    GAMMA ( JlfNS)+BET A(J , NS1J/C 

KK=NS/2 
C 
C 

C         WRITE  THE  CURRENT  INDUCED  AT  THE  CENTER  OF  THE  WIRE 
C 
C 

WRITE16, 96901  J,BETA(J,KK) 
9690   FORMATUOXf  I4f21Xf  E12.5) 
200    CONTINUE 
C 

C  END  OF  CURRENTS  COMPUTATION 

C 
C 
C 

C         START  OF  SCATTERED  FIELD  COMPUTATION 
C 
C 

DO  400  LP=1,ITS 

ES(LP)=0.0 

DO  500  K=1,NS1 
IF  (TS.GT.90.0)  GO  TO  501 

IF  (TS.EQ.90.0)  GO  TO  503 
C 

C         FOR  REFRACTION  ANGLE  LESS  THAN  90  DEGREES 
C 

KP  =  K 

NQ=LP-1-(KP-0.5I*AR 

NQP=LP-l-(KP+0.5)*AR 

GO  TO  502 
C 

c 
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C         FOR  REFRACTION  ANGLE  MORE  THAN  90  DEGREES 
C 

501  CONTINUE 
KP=K-NS 

NQP=LP-1- (KP-0.5 ) *AR 

NQ=LP-L-(KP+0.5)*AR 

GO    TO    502 
C 

C  FOR    REFRACTION    ANGLE    EQ.    90    DEGREES. 

C 
503         CONTINUE 

NQ=LP-1 

NQP=NQ 

502  CONTINUE 
NQ1=NQ-1 
NQ2  =    NQP-1 

IF(NQ.LE.O)    GO   TO    500 
IFINQ.EQ. 1)    GO   TO    510 
IF(NQ.EQ.NQP)    GO    TO    520 
ZQ=(LP-1-NQ)/AR    +0.5 
EPS=DZ*(ZQ-KP) 
DEL=DZ-EPS 

ES(LP)=ESILP)+EPS*<BETA(NQfK)-BETA(NQl,K) J 
IFCNQ.EQ.2)    GO   TO    530 

ES(LP)=ES(LP) +DEL*     <BETA< NQP, K)-BETA(NQ2, K ) ) 

GO    TO    500 
510         CONTINUE 

ES(LP)=ES(LP)+DZ*2*    BETA(NQ,K) 

GO    TO    500 
520         CONTINUE 

ES{LP)=ES(LP)+DZ*(BETA(NQ,KJ-BETA(NQl,K)) 

GO    TO    500 
530         CONTINUE 

ES(LP)=ES(LP)+DZ*2*    BETA<NQP,K) 
500         CONTINUE 

CT=-SIN(TSR)*<1.0    E-71/DT 

ES(LP)=CT*ES<LP) 
400         CONTINUE 
C 
C 
C 

C  END    OF    FIELD    COMPUTATION 

C 
C 

c 

c 

C         PLOT  OF  CURRENT  AT  CENTER  OF  THE  WIRE 

C  PLOT  EVERY  MM  POINT. 

C 

MM=ITS/iOO 

00  600  I=MM,ITS,MM 

K=I/MM 

X(K1  =I*DZ 

KK=NS/2 

Y(K)  =  BETA(ItKK)  *1000.0 
600    CONTINUE 

MM=ITS/80 

WRITE(&,9630) 
9630       FQRMAT(lHl) 

MODCUR=0 

CALL    PL0TP<X,Y,80    ,MODCURI 

WRITE    (6,9600)       TI 
9600       FORMAT*// //,5X,« CURRENT    INDUCED    AT    CENTER    OF    WIRE1,' 

2  SCATT£RER«,    /,5X, 'FIELD   INCIDENT    AT    ANGLE' t lF6.1t///t 

3  lOXt'X-AXlS    TIME     IN    LIGHT-METERS1, 

4  //,10X,«Y-AXIS    CURRENT    IN    MILLIAMPSM 
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c 

C  PLOT    OF    SCATTERED    FIELD 

C 

DO    700    I=MM,ITS,MM 

K=I/MM 

Y(K)=ES(I)  *1000.0 
700    CONTINUE 

WRITE  (6,9640) 
9640   FORMATUHi) 

MODCUR=0 

CALL  PLOTP(X,Y,80  ,MODCUR) 

WRITE    (6,9620)    TS 
9620      F0RMAT(////,5X, 'NORMALIZED    FAR    FIELD    8ACK    SCATTERED', • 
2BY    WIRES     /,5X,'AT    ANGLE* , 1F6 .It / / f 10X, • X-AX I S    TIME1,' 
3       LIGHT    METERS' ,//,10X, 'Y-AXIS    FAR    FIELD    MILLIVOLTS') 
C 
C 

C  END    OF    PRCGRAM    GO    BACK    TO    NEXT    ANGLE    IF    ANY 

C 

GO    TO    10 
999         CONTINUE 

STOP 

END 
C 
C 
C 

c 

c 
c 

C         EXAMPLE  OF  DATA  CARDS 

C 

C 

C 

90.0   90.0 
1.0 
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APPENDIX  E 
C 
C 
C 

c 
c 

c 
c 
c 

C      ************************************************ 

c  *  * 

C  *  TIME  DOMAIN  RESPONSE  AND  INVERSE  SCATTERING  * 

C  *  * 

C  *    FOR    A    PERFECT    CONDUCTOR    BODY    OF    REVOLUTION  * 

C  *  * 

c  ************************************************ 

c 
c 
c 

c 
c 

C  NAME    OF    PROGRAM    :    GAIL 

Q  **************** 

C 

c 
c 

c 
c 

C  PURPOSE: 

Q  ******** 

C 

C  THIS    PROGRAM    COMPUTES    THE    TRANSIENT    TIME    DOMAIN 

C  RESPONSE    OF    A    PERFECT    CONDUCTOR    BODY   OF    REVOLUTION* 

C  AND    APPLIES    IT    TO    THE    SOLUTION    OF    THE    INVERSE 

C  SCATTERING    PROBLEM. <CCMPUTATI ON    OF    THE    SHAPE    OF    THE 

C  BODY    FRO^I    ITS    TIME    DOMAIN    RAMP    RESPONSE) 

C 

C 

c 
c 
c 
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C  USAGE: 

Q  ****** 

C 

C  THE  PROGRAM  CAN  BE  USED  IN  ONE  OF  FOUR  WAYS 

C  DEPENDING  OF  THE  DESIRED  OUTPUT  AND  CONSEQUENTLY  OF 

C  THE  VALUE  ASSIGNED  TO  THE  INITIATING  PARAMETERS  'INC* 

C  AND  « INVER». 
C 

C  THE  POSSIBLE  USES  ARE: 
C 

C  1.  COMPUTATION  OF  THE  TRANSIENT  TIME  DOMAIN  IMPULSE 

C  RESPONSE  (BACKSCATTERED  FIELD  ),WHEN  THE  INCIDENT 

C  FIELD  IS  A  GAUSSIAN  IMPULSE  .(INC=1) 

C  THE  CURRENTS  INDUCED  ON  THE  BODY  CAN  BE  MADE  AVAILABLE 

C  AT  THE  OUTPUT  IF  DESIRED  (FOR  EMP  PROBLEMS). 
C 

C  THE  INPUTS  REQUIRED  ARE: 
C 

C  A.  THE  SPACE  SAMPLES  OF  THE  CONTOUR  FUNCTION  OF 
THE  BODY,RO(Z),ENTERD  IN  DATA  CARDS. 


I 


C  B.    THE    NUMBER    OF    SEGMENTS    SUBDIVIDING    THE    CONTOUR 

C  OF    THE    BODY    (NS<20) . 

C 

C  C.  THE  LENGTH  OF  A  TIME  INCREMENT  (DT). 

C 

C  D.  THE  NUMBER  OF  TIME  SAMPLES  NEEDED  ( ITS<=  60). 

C 

C  E.  THE  LENGTH  OF   EACH  SEGMENT  DS(20). 

C 

C  F.  THE  PARAMETERS  DEFINING  THE  GAUSSIAN  IMPULSE. 

C 

C        2.  COMPUTATION  OF  THE  RAMP  RESPONSE  OF  THE  80DY 

C      WHEN  THE  INCIDENT  FIELD  IS  A  RAMP.(INC=2  AND  INVER= 1 ) 

C 

C      THE  INPUTS  REQUIRED  ARE: 

C 

C  SAME  AS  A  THRU  E  LISTED  ABOVE. 

C 

C        3.  INVERSE  SCATTERING  SOLUTION  STARTING  FROM  THE 

C      BODY  ITSELF.  (THIS  CASE  IS  A  SIMULATION  PROCEOURE, 

C      WHERE  THE  PROGRAM  RETRIEVES  BY  INVERSE  SCATTERING 

C      THE  SHAPE  OF  THE  KNOWN   BODY.  (INC=2  AND  INVER=1). 

C 

C      THE  INPUTS  REQUIRED  ARE: 

C 

C  SAME  AS  A  THRU  E  LISTED  ABOVEtPLUS  THE  LENGTH  OF 

C      A  SEGMENT  FOR  NEW  SUBDIVISION  OF  THE  BODY.(DSO) 

C 

C        4-.  INVERSE  SCATTERING  SOLUTION  STARTING  FROM  THE 

C      MEASURED  GIVEN  RAMP  RESPNSE  OF  AN  UNKNOWN  BODY  (HSFM). 

C      THIS  CASE  IS  THE  TRUE  COMPUTATION  OF  THE  SHAPE  OF  THE 

C      BODY  WHEN  ONLY  ITS  RAMP   RESPONSE  IS  KNOWN. 

C      (INC=2  AND  INVER=2) 

C 

C      THE  INPUTS  REQUIRED  ARE: 

C 

C  A.    THE    TIME    SAMPLES    OF    THE    RAMP    RESPONSE . (ENTERED 

C  IN    DATA    CARDS) 

C 

C  B.  THE  LENGTH  OF  THE  DESIRED  SEGMENTS.  ( D SO ) 

C 

C  C.  B  THRU  D  LISTED  ABOVE. 

C 

C 

C        NOTE:  NO  INVERSE  SCATTERING  IS  COMPUTED  WHEN  INC=l  . 

C 
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C  PARAMETER    DEFINITION: 

c 

C      INC:  DEFINES  THE  TYPE  OF  INCIDENT  FIELD 

C 

C  INC=1    THE    INCIDENT    FIELD    IS    A    GAUSSIAN    IMPULSE. 

C  INC=2    THE     INCIDENT    FIELD    IS    A    RAMP. 

C 

c 

C  INVER  :   DEFINES  THE  TYPE  OF  PROGRAM  REQUIRED 

C 

C  INVER=0    NO    INVERSE    SCATTERING 

C  INVER=l       THE    PROGRAM   COMPUTES    THE     INVERSE    SCATTERING, 

C  STARTING    FROM    A    GIVEN    BODY    SHAPE. 

C  INVER=2       THE    PROGRAM    COMPUTES    THE     INVERSE     SCATTERING 

C  FROM    A    GIVEN    RAMP    RESPONSE    OF    THE    BODY. 

C 

C 

C  NS:  NUMBER  OF  SEGMENTS  SUBDIVIDING  THE  CONTOUR  (<=20) 

C 

C  ITS:  NUMBER  OF  TIME  SAMPLES  NEEDED  ( <=  60) 

C 

C  DT:  LENGTH  OF  THE  TIME  INCREMENT  (SECONDS) 

C 

C  DS(20)  :  LENGTH  OF  THE  SEGMENTS  SUBDIVIDING  THE 

C  CONTOUR  OF  THE  BODY. ( MAXIMUM  20  SEGMENTS) 

C  TO  BE  ENTERD  IN  DATA  STATEMENT. 

C 

C  2(21):  Z-COORDINATE  OF  THE  SPACE  SAMPLE  POINTS  ON  THE 

C  CONTOUR  OF  THE  BODY. ( MAXIMUM  21  POINTS  ENTERED 

C  DELIMITING  20  SEGMENTS ).( METER) 

C 

C  RO  (21)  :RADIUS  OF  THE  BODY  AT  THE  SPACE  SAMPLE  POINTS 

C  DEFINED  BY  Z.(MAX.  21  ). (METER) 

C 

C  CJS(20,60):  SPACE  TIME  SAMPLES  OF  THE  S  COMPONENT  OF 

C  THE  CURRENTS  INDUCED  ON  THE  BODY.  (UP  TO 

C  20  SPACE  SAMPLES  AND  60  TIME  SAMPLES  EACH) 

C 

C  CJPHI(20,60):  SPACE  TIME  SAMPLES  OF  THE  PHI  COMPONENT 

C  OF  THE  CURRENTS  INDUCED  ON  THE  BODY.  (UP  TO 

C  20  SPACE  SAMPLES  AND  60  TIME  SAMPLES  EACH) 

C 

C  HSF  (60):   TIME  SAMPLES  OF  THE  SCATTERED  FIELD. UP  TO 

C  60  SAMPLES. 

C 

C  HSFM(60):       TIME    SAMPLES    OF    THE    INITIAL    RAMP    RESPONSE. 

C  (     THE    MEASURED    RAMP    RESPONSE) 

C 

C  HSFC(60):  TIME  SAMPLES  OF  THE  COMPUTED  RAMP  RESPONSE 

C  OF  THE  BODIES  OBTAINED  BY  INVERSE  SCATTERING 

C 

C  NP(20):  NUMBER  OF  PATCHES  ON  THE  RINGS  DEFINED  BY 

C  EACH  SEGMENT. 

C 

C  B(20):  ANGLE  SUSTAINED  BY  THE  PATCHES  ON  EACH  RING. 

C 

C  ALFAK:  SLOPE  OF  THE  CONTOUR  OF  THE  BODY  CONTOUR 

C  FUNCTION  tMEASURED  AT  THE  CENTER  OF  THE 

C  EACH  SEGMENT. 

C 

C  SI     (20):    SIN(ALFAK)     FOR    EACH    SEGMENT. 

C 

C  CO  (20):  COS(ALFAK)  FOR  EACH  SEGMENT. 

C 

C  EPS:  NORMALIZED  MEAN  SQUARED  ERROR  ON  THE  BODY  SHAPE 

C  COMPUTED  BY  INVERSE  SCATTERING  PROCEDURE. 
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c 

C      ITER:  NUMBER  OF  ITERATIONS  EXECUTED  FOR  INVERSE 

C  SCATTERING  SOLUTION. 

C 

C      C:  SPEED  OF  LIGHT  (METER/SEC) 

C 

C      CP  :  CONSTANT  C*DT. 

C 

C      OTHER  PARAMETERS  USED  ARE  DEFINED  IN  THE  PROPER 

C      SUBROUTINE. 

C 

C 

C 

C        REMARKS: 

£  ******** 

C 

C  1.    THE    SHAPE    OF    THE    BODY    IS    ENTERD    IN    DATA 

C  CARDS    WITH    THE    FORMAT    2F10.5 

C 

C  2.  THE  TIME  SAMPLES  OF  THE  MEASURED  RAMP 

C  RESPONSE  HSFM  ARE  ENTERD  IN  DATA  CARDS  WITH  THE 

C  FORMAT   1F10.5 

C  THE  VALUES  OF  THE  SAMPLES  MUST  BE  NORMALIZED  THAT 

C  MEANS  MULTIPLIED  BY  C  AND  R  (RANGE  TO  TARGET) 

C 

C  3.  IN  THIS  PROGRAM  THE  INCIDENT  FIELD  H  IS 

C  ASSUMED  TO  BE  VERTICAL  (HORIZONTAL  POLARIZATION). 

C  THE  DIFFERENCE  BETWEEN  THE  RESPONSE  TO  THE  TWO 

C  POLARIZATIONS  IS  ONLY  THE  SIGN  OF  THE  RESPONSE. 

C  (  SEE  ALSO  REMARKS  IN  SUBROUTINE  SHAPE) 

C 

C 

c 

C        SUBROUTINES  USED: 
C        ***************** 

c 

C  1.  SUBROUTINE  •FIELD1  :  COMPUTES  THE  BACK  SCATTERED 

C  FIELD  FOR  BOTH  RAMP  AND  IMPULSE  EXCITATION. 

C 

C  2.  SUBROUTINE  'SHAPE'  :  COMPUTES  WITH  A  METHOD  OF 

C  ITERATIONS,  THE  SHAPE  OF  THE  BODY  FROM  ITS  RAMP 

C  RESPONSE. (SEE  SUBROUTINE  SHAPE  FOR  DETAILS.) 

u 

C 

C        3.  LIBRARY  SUBROUTINE  'PLOTP'  FOR  PLOTS  OF  FIELDS 

C        AND  SHAPE  OF  THE  BODY. 

C 

C 

C 

C        METHOD  OF  SOLUTION: 

Q  ******************* 

C 

C  SEE  SUBROUTINE  FIELD  FOR  THE  COMPUTATION  OF  IMPULSE 

C  AND  RAMP  RESPONSES. 

C 

C  SEE  SUBROUTINE  SHAPE  FOR  THE  COMPUTATION  OF  THE 

C  SHAPE  OF  THE  BODY. 

C 

C  INVERSE  SCATTERING:  THE  FIRST  APPROXIMATION  OF  THE 

C  SHAPE  IS  OBTAINED  FROM  ITS  MEASURED  RAMP  RESPONSE. 

C  A  NEW  RAMP  RESPONSE  IS  COMPUTED  FOR  THAT  APPROXIMATE 

C  BODY. THE  PHYSICAL  OPTICS  PART  OF  THE  RESPONSE  IS 

C  EXTRACTED, AND  IS  USED  TO  COMPUTE  A  NEW  SHAPE, AND  SO 

C  ON  ,UNTILL  A  MINIMUM  MEAN  SQUARED  ERROR  IS  OBTAINED 

C  BETWEEN  THE  LAST  COMPUTED  RESPONSE  AND  THE  MEASURED. 

C  (SEE  ALSO  SECTION  5  OF  THIS  THESIS) 

C 
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c  *  * 

C  *    MAIN    PROGRAM       * 

C  *  * 

C 

C  SPHERE  EXAMPLE 


C 


DIMENSION  X(60) 

DIMENSION  CJS<20,60   I ,CJPHI ( 20 ,60) ,  HSF(60) 


DIMENSION  Z(21),R0(21) ,NP(20) 

DIMENSION  SI  (20)  ,C0  ( 20 ) » B  (20)  »DS  (20) 

DIMENSION  HSFM(60)  ,HSFC  (60) 

COMMON    HSF,Z,R0,ITS,NS,NS1 

COMMON    XMAX 

DATA    HSFM/60*0.0/ 

DATA    HSFC/60*0.0/ 

DATA    DS/20*0. 20900/ 
C 
C 

C        START  PROGRAM 
C 

C        DEFINE  TYPE  OF  PROGRAM, INVER  AND  INC 
C 

INC=2 

INVER=1 
C 

C        COMPUTATION  PARAMETERS 
C 

NS  =  15 

NS1=NS+1 

ITS=60 

IT=ITS 

DSO=0.209 

C=3.0E+8 

DT=0.2/C 

C?=C*DT 
C 

IF  (  INC.EQ.l)  GO  TO  LI 

IF  (INVER   .EQ.2)  GO  TO  15 
C 

C        STARTING  POINT   FOR  INV£R=1  AND  INVER=0 
C 

11     CONTINUE 
C 

C        READ  THE  SHAPE  OF  THE  BODY 
C 

DO    10    K=1,NS1 

READ(5,9510)     Z(K) ,RO(K) 
9510       F0RMAT(2F10.5) 
10  CONTINUE 

C 

C  PLOT    THE    INITIAL    SHAPE    OF    THE    BODY 

C 

NS2=NS1+1 

XMAX      =    Z (NS1)*1.25 

Z(NS2)=    XMAX 

RO    (NS2)=    Z(NS2) 

WRITE(6,9620) 

9620  FORMAT    UH1) 
MODCUR=0 

CALL    PLOTP< Z,R0,NS2,M0DCUR) 
WRITE    (6,9621) 

9621  FORMAT    { /// , 20X , ' I  NIT  I AL    SHAPE    OF    BODY' , // ,21X , 
2       • X-AXIS    Z,    Y-AXIS    RADIUS     (R0)«) 

C 
C 
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C  COMPUTE   THE    BACK-SCATTERED    FIELD 

C 

C 

CALL    FIELD    (DS»CP , DT, INC ) 

c 

C  END    OF    PROGRAM    FOR    INVER=0    AND    INVER=1 

C 

IF<  INC.EQ.l)     GO    TO       99 

IF     (INVER       .EQ.O)     GO    TO   99 

IF    (INVER       .EQ.il    GO   TO    L9 
15  CONTINUE 

C 

C  STARTING    POINT    FOR    INVER=2 

C 

IN02 
C 

C        READ  THE  GIVEN  RAMP  RESPONSE 
C 

DO  12  N=i,ITS 

READ  (5,9520)  HSFM  (N) 
9520   FORMAT  (1F10.5) 

HSF(N)=HSFM(N) 
12     CONTINUE 

GO  TO  920 
C 

19     CONTINUE 
C 

C        STORE  THE  COMPUTED  RAMP  RESPONSE 
C 

DO  910  N=l,ITS 

HSFM(N)=HSF(N) 
910    CONTINUE 
920    CONTINUE 
C 

C        STARTING  POINT  FOR  INVERSE  SCATTERING  COMPUTATION 
C 

ITER=1 
C 

C        COMPUTATION  OF  THE  FIRST  APPROXIMATION  OF  THE  SHAPE 
C 

c 
c 

CALL   SHAPE  { OS, ITER,CP ,DSO, I TJ 
C 

C 

c 

C        COMPUTE  SUM  OF  SQUARES  OF  THE  TIME  SAMPLES  OF  THE 

C        ORIGINAL  RAMP  RESPONSE. 

C 

ITS=IT+10 

SIGMA1=0.0 

00  940  N=l,  ITS 

SIGMA1=SIGMA1+HSFM(N)**2 
940     CONTINUE 
C 

950     CONTINUE 
C 

C        COMPUTE  THE  RESPONSE  OF  THE  APPROXIMATE  BODY. 
C 
C 
C 


C 

C 

c 
c 


CALL    FIELD    ( DS,CP , DT, INC ) 
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C        COMPUTE  NORMALIZED  MEAN  SQUARED  ERROR. 
C 

SIGMA2=0.0 

0-3  960  N=lf  ITS 

SIGMA2=SIGMA2+(HSFM(N)-HSF(N) )**2 

960  CONTINUE 
EPS=SIGMA2/SIGMA1 

C 

C        WRITE  ERROR 

C 

WRITE  (6,9666)  EPS 
9666    FORMAT  {   //, 15X ,  » NORMALIZED  MEAN  SQ.  ERROR8  ,1F10. 5) 

IF  (ITER.GT.l)  GO  TO  961 

EPS1=1. 

961  CONTINUE 
C 

C        END  THE  PROGRAM  WHEN  THE  ERROR  IS  MINIMUM 
C 

IF(EPS.GT.EPSl)  30  TO  999 

EPS1=EPS 
C 

C        CONTINUE   WHEN  THE  ERROR  IS  NOT  YET  MINIMIZED 
C 

ITER=ITER+1 
C 

C        EXTRACTION  OF  THE  PHYSICAL  OPTICS  RESPONSE 
C 

WRITE  (6,9649) 
9649   FORMAT  (1H1,  10X, • PHYSI CAL  OPTICS  RESPONSE »,//, 9X ,» N« , 
2   10X,«H(N) •,//) 

DO  970  N=1,ITS 

HSFC(N)=HSFC(N)+HSF(N) 

HSF(N)=ITER*HSFM(N)-HSFC(N) 
WRITE  (6,9651)  N,HSF(N) 
9651   F0RMAT(5X,I5,5X,E12.5) 
970     CONTINUE 
C 

C        USE  PHYSICAL  OPTICS  RESPONSE  TO  COMPUTE  A  NEW  SHAPE 
C 
C 

CALL   SHAPE  (  DS  ♦  I  TER,CP  ,DSO,  I  T  ) 
C 

r 

C        RETURN  TO  COMPUTE  THE  RESPONSE  TO  THE  NEW  30DY 
C 

GO  TO  950 
999  CONTINUE 
C 

C        END  OF  PROGRAM  WHEN  MINIMUM  ERROR  ACQUIRED 
C        WRITE  FINAL  ERROR 
C 

WRITE  (6,9625)  ITER.EPS1 
9625   FORMAT  ( // ,20X, • F INAL  SHAPE  AFTER ',  14, ' ITERATIONS ■ ,/// 

2   ,20X, 'NORMALIZED  MEAN  SQUARED  ERROR  IS»,1F10.5) 
99       CONTINUE 
C 

C        END  OF  PROGRAM 
C 

STOP 

END 
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c 
c 
c 

c 

c      ************************ 
C      *  * 

C      *    SUBROUTINE  FIELD    * 
C      *  * 

C 

c 
c 
c 

C        PURPOSE: 
C        ******** 

C 
C 

c 

C  THIS  SUBROUTINE  PROVIDES  THE  SOLUTION  TO  THE 

C  MAGNETIC  FIELD  INTEGRAL  EQUATION  APPLIED  TO  A 

C  BODY  OF  REVOLUTIONtEXCITED  BY  A  GAUSSIAN  PULSE 

C  OR  A  RAMP  FIELDS. THE  CURRENTS  INDUCED  3N  THE  BODY 

C  AND  THE  BACK-SCATTERED  FIELD  ARE  COMPUTED. 

C 

C  USAGE: 

Q  ****** 

c 
c 

C        CALL  FIELD  (DS,  CP  t  DT,  INC) 

C 

c 

C        PARAMETER  DEFINITION: 
-        ********************* 

C 

c 

c 

C        HI:  INCIDENT  FIELD 

C 

C        C3,C4:  CONSTANTS. 

C 

C 

C        METHOD  OF  SOLUTION: 

Q  ******************* 

c 

c 

C  THE    METHOD    OF    SOLUTION    IS    DESCRIBED    IN    SECTION    4 

C  OF   THIS    THESIS. 

C 

C 

c 
c 
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c 
c 


SJBROUTIME    FIELD    ( DS, CP »DT, INC ) 


DIMENSION  X(60) 

DIMENSION  CJS(20T60   ) , CJPHI ( 20 , 60 ) ,  HSF(60) 

DIMENSION  Z(21),RQ<21) ,NP(20) 

DIMENSION    SI     (20)     ,  CO    ( 20  ) ,  B    (20)     ,DS    (20) 

COMMON    HSF,ZtRO,ITS,NSTNSl 

c 

C        CONSTANT  DEFINITION 
C 

PI=3. 141592654 

C=3.0E+8 

C4=-1./(DT*4.0*:> 

C3=l./( PI*2.0) 
C 
C  GAUSSIAN    IMPULSE    PARAMETERS 


C 


c 


AA=0.6    E+9 
TMAX=DT*5.0 
A2  =  AA*AA 


DO    20    K=1,NS 

K1=K+1 
C 

C  COMPUTE    THE    SLOPE    OF    THE    SEGMENTS 

C 

WK=(R0(K1 )-RO(K))/(Z(Kl)-Z(K) ) 

ALFAK=ATAN( WK ) 

SI(K)=SIN(ALFAK) 
CC( K)=COS(ALFAK) 
C 

C  COMPUTE    THE    CENTRAL    SPACE    SAMPLE    POINTS 

r 

Z(K)=(Z(K1)+Z(K) J/2.0 

RO < K) = ( RO ( K 1 ) +R0 ( K ) ) / 2 . 0 

C 

C        COMPUTE  THE  NJM3ER  OF  PATCHES  IN  EACH  RING 

C 

PK=     PI*RO(K)/DS(K) 

NP(K)=2*IFIX(PK) 

IF  (NP(K) .GT.O)  GO  TO  23 

NP(K)=i 
C 

C        COMPUTE  THE  ANGLE  SUSTAINED  BY  THE  PATCHES  OF  EACH 
C        RING. 
C 

23  CONTINUE 
B(K)=2.0*PI/NP(K) 

C 

C        INITIAL  CONDITIONS  OF  THE  CURRENTS 

TZK=1.0+Z(K)/CP 

DO  25  N=l,2 
IF  (N.LT.TZK)  GO  TO  26 
IF(  INC.EQ.l )  GO  TO   24 
C 

C        WHEN  INCIDENT  FIELD  IS  A  RAMP 
C 

HI=(N-1)  *CP-Z(K) 
GO  TO  27 

24  CONTINUE 
C 
C        WHEN  INCIDENT  FIELD  IS  AN  IMPULSE 


ARG=(N-I) *DT~TMAX-Z(K)/C 
HI=(N-1)*EXP(-A2*ARG*ARG) 
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GO  TO  27 

26  CONTINUE 
HI=0.0 

27  CONTINUE 
C 

C  INITIAL    CURRENTS    INDUCED 

C 

CJPHI(K,N)=-2.0*HI*SI(K) 

CJS(K,N)=-2.0*HI 
25       CONTINUE 
20     CONTINUE 
C 

c 

C        COMPUTE  THE  CURRENT   COMPONENTS 

C 

C        COMPUTE  FIRST  THE  INTEGRAL  PART   OF  THE  EQUATION 

C        AT  TIME  N 

C 

DO  100  N=3,ITS 
C 

C        AT  THE  K-TH  SAMPLE 
C 

DO  200  K=l,NS 

CJS(K,N)=0.0 

CJPHI(K,N)=0.0 
C 
C        INFLUENCE  OF  ALL  OTHER  SPACE  CURRENTS 


C 


c 


c 
c 


DO  300  1=1, NS 
NP2=NP(I) 


DO  400  M=l,NP2 

IF  (M.EQ.l)  GO  TO  4-10 
GO  TO  420 
410  CONTINUE 

IF  (I .EQ.K)  GO  TO  400 
420         CONTINUE 
C 

C        COMPUTE  DISTANCE  BETWEEN  SAMPLES 
C 

CM=COS( (M-1)*B(I) ) 

SM=SIN( (M-1)*B(I) ) 
V=Z(I)-Z(K) 

R2=(R0(K) *ROt  Kl+ROC I ) *R0( I ) +V*V-2. 0*R0< I > *R0 ( KJ *CM) 
R=SQRT(R2) 
C 

C  COMPUTE    RETARDATION    TIME 

C 

Q=N-R/CP 
C 

C  INTERPOLATION    OF    THE    CURRENT    AND    THEIR    TIME 

C  DERIVATIVE    AT    THE    RETARDED    TIME 


NQ=IFIX(Q) 

NQP=NQ+1 

NQM=NQ-1 

NQ2=NQ+2 

T1=Q-NQ2 

T2=Q-NQP 

T3=Q-NQ 

T4=Q-NQM 

IF     (Q,LE.1.0)    GO    TO    400 

IF    (Q.LE.2.0)    GO    TO    435 

IF    (N.EQ.NQP)    GO    TO    430 

IF     (N.EQ.NQ2)    GO    TO    432 
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C  FCUR    POINT    LAGRANGIAN    INTERPOLATION 

C 

C 

Yl=T2*T3*T4*CJS(ItNQ2)/6.0-T3*T4*Ti*CJS< I, NQP)/ 2.0 
2         +T4*T1*T2*CJS(  I»NQ)/2.0-Tl*T2*T3*CJS(  I ,NOM)/6.0 

Y2=T2*T3*T4*CJPHI( 1 1 NQ2  )/6. 0-T3*T4*Tl*C JPHI < I, NQP)/ 2.0 
2         +T4*Ti*T2*CJPHI( I,NQ)/2.0-Tl*T2*T3*CJPHIU ,NQM)/6.0 
DJS=(T2*T3+T2*T4+T3*T4)*CJS<  I,NQ2)/6.0 

2  -(T3*T4+T4*TI+T1*T3)*CJS<  I»NQP)/2.0 

3  *(T4*TH-T4*T2+TL*T2)*CJS(  ItNQ)/2.0 

4  -(T1*T3+T2*T3+T2*T1)*CJS< I»NQM)/6.0 
DJPHI=(T2*T3+T2*T4+T3*T4)*CJPHI( I, NQ2I/6.0 

2  -<T3*T4+T4*T1+T1*T3)*CJPHI(  I,NQP)/2.0 

3  +<T4*T1+T4*T2+T1*T2)*CJPHI( I,NQ)/2.0 

4  -(T1*T3*T2*T3+T2*T1)*CJPHI(I fNQMI/6.0 
GO    TO    440 

430         CONTINUE 

C 

C 

C        TWO    POINT  LAGRANGIAN  INTERPOLATION 

C 

c 

Y1  =  T4*CJSU,NQ    )-T3*CJS(IfNQMJ 
Y2=T4*CJPHI(I TNQ    )-T3*CJPHI(I ,NQMI 
DJS=CJS(I,NQ    )-CjS(I?NQMI 
DJPHI=CJPHI(I  ,NQ    l-CJPHKI,  NQN1 
GO    TO    440 

432  CONTINUE 
C 

C 

C        THREE  POINT  LAGRANGIAN  INTERPOLATION 

C 

C 

Y1=T4*T3*CJS( I,^QP»/2.0-T2*T4*CJS(ItNQ    ) 
2      +T2*T3*CJS( I.NQMI/2.0 

Y2=T4*T3*CJPHIU,MQP)/2.0-T2*T4*CJPHI(I,NQ    ) 
2      +T2*T3*CJPHl(I,NQM)/2.0 

0JS=(T4+T3)*CJS(I »NQP»/2.0-(T2+T4)*CJS(I,NQ) 
2  +  <T2*T3)*CJS(  I»NQM)/2.0 

DJPHI=<T4+T3)*CJPHI(I,NQP)/2.0-{T2+T4)*CJPHI( IfNQ    ) 
2  ♦{T2*T3)*CJPHI(I,NQM»/2.0 

GO    TO    440 
435         CONTINUE 

IF    (N.EQ.NQP)    GO    TO    400 
IF    (N.EQ.NQ2)    GO    TO    433 
C 

C  THREE    POINT    LAGRANGIAN    INTERPOLATION 

C 

Y1=T2*T3*CJS( I,NQ2)/2.0-T3*Tl*CJS( I f NQP I 
2      +T2*T1*CJS( I,NQ)/2.0 

Y2=T2*T3*CJPHI(I,NQ2)/2.0-T3*Tl*CJPHI(I,NQP) 
2      +T2*Tl*CJPHI(I,NQ)/2.0 

0JS=(T2+T3)*CJS( I,NQ2)/2.0-(T3*Tl)*CJS< I, NQP) 
2  +(T2+T1)*CJS( ItNQI/2.0 

DJPHI=(T2+T3)*CJPHI<I,NQ2)/2.0-<T3+Tl)*CJPHI( I, NQP) 
2       +(T2+Tl )*CJPHI(I,NQ)/2.0 
GO  TO  440 

433  CONTINUE 
C 

C        TWO    POINT  LAGRANGIAN  INTERPOLATION 
C 

Y1=T3*CJS( I,NQP)-T2*CJS(I,NQ) 

Y2=T3*CJPHI(I ,NQP)-T2*CJPHI (I ,NQ) 

DJS=CJS(I,NQP)-CJS(I,NQ) 

DJPHI=CJPHI(I,NQP)-CJPHI(I,NQ) 
C 
440         CONTINUE 
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c 

C        FINAL  COMPUTATION  OF  THE  INTEGRAL 
C 

v=zm-z(Kj 

F1=V*CM*SI<  n*(R0(K)-R3(I)*CM)*C0U) 

Xl=Fl*CM 
X2=V*(SM*S^) 
U=V*SI(  I}*SI(K)+RO(K)*CO<K)*SI( I )-RO< I ) *SI < Kl *CO< I ) 
F2=(R0(K)*CM-R0( I) )*CO( K) +CM*V*S  I  < K) 
Y11=Y1*XI 
Y22=Y2*X2 
XPl=F2*CM 
XP2=U*(SM*SM) 
YP1=Y2*XP1 
YP2=Y1*XP2 
C1  =  DS( I)*RO(  I)*B(I  )/(CP*R2) 
R3=R*R2 

C2=OSl  I)*  RO( I )*B(I)/R3 
CJPHI(K,N)=CJPHI<K,N)+C1*(DJPHI*XP1+DJS*XP2) 
2    +C2*< YP1+YP2) 
CJS(K,N)=  CJS<K,N)+C1*(DJPHI*X2  >DJS*X1  ) +C2* < Y11+Y22 ) 
400         CONTINUE 
C 
300        CONTINUE 

CJS(K,N)=CJS(K,NI*C3 
CJPHU K,N)=CJPHI(K,NI  *C3 
C 

c 

C        ENO  OF  COMPUTATION  OF  THE  INTEGRAL 
C 

c 

C        ADDITION  OF  THE  CURRENT  INDUCED  8Y  THE  INCIDENT 

C        FIELD  AT  THE  COMPUTATION  POINT . (PHYSICAL  OPTICS  PART 

C 

TZK=1.0+Z(K)/CP 

TD=TZK+2.*TMAX/DT 

IF  (N.LT.TZK)  GO  TO   200 

IF( INC.EQ.l)  GO  TO   240 
C 

C        WHEN  INCIDENT  FIELD  IS  A  RAMP 
C 

HI  =  (N-1)*CP-ZU) 

GO  TO  250 
240    CONTINUE 

IF  (N.GT.TD)  GO  TO  200 
C 

C  WHEN    INCIDENT    FIELD    IS   AN    IMPULSE 

C 

ARG=(N-1) *DT-TMAX-Z(K)/C 
HI=EXP(-A2*ARG*ARG) 
250         CONTINUE 
C 

C  FINAL    CURRENT    AT    THE    K-TH    SPACE    SAMPLE    AT    TIME    N 

C 

CJS(K,N)=CJS(K,N)-2.0*HI 
CJPHI{K,N)=CJPHI(K,N)-2.0*HI*SI{K1 
C 

C  NEXT    SPACE    SAMPLE 

C 

200  CONTINUE 

C 

C  NEXT    TIME 

C 

100    CONTINUE 
C 

C        END  OF  COMPUTATION  OF  ALL  SPACE-TIME  CURRENTS 
C 
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IF    (INC.EQ.l)    GO    TO    110 

WRITE(6t9623) 
9623      FORMAT( lrl  1 ,  10X,  •  NORMAL  I  ZED    FIELD       (H*C*R)'t//» 
2       10X,     'RAMP    RESPONSE1  ,//,9X,  •N'  ,10X,»H(M)«  ,//) 
GO    TO    111 

110  CONTINUE 
WRITEI6, 97231 

9723       FORMAT! 1H1,10X, 'NORMALIZED    FIELD       (H*R)        •,//, 

2       10X,    'IMPULSE    RESPONSE* »//,9X,»N» v10XvaH{N) ' »//) 
C 
C 

C  START    COMPUTATION    OF    THE    BACK-SCATTERED    FIELD 

C 
C 

111  CONTINUE 
C 

C        AT  THE  N-TH  TIME 
G 

DO  500  N=l,ITS 
HSF(N)=0.0 

IF  (N.EQ. 1)  GO  TO  650 
C 

C        FROM  EACH  SAMPLE 
C 

DO  600  K=1,NS 
C 

C        COMPUTATION  OF  THE  RETARDATION  TIME 
C 

Q=N-Z(K)/CP 
NQ=IFIX(Q) 
C 

C  INTERPOLATION    OF    THE    TIME    DERIVATIVE    OF    THE    CURRENT 

C  AT    THE    RETAROED    TIME 

C 

NQP=NQ+1 

NQM=NQ-1 
NQ2=NQ+2 
T1=Q-NQ2 
T2=G-NQP 
T3  =  Q-NQ 
T4=Q-NQM 

IF    (Q.LE.1.0)    GO    TO    600 
IF    <Q.LE.2.0)    GO    TO    601 

IF     (N.EQ.NQP)    GO    TO    602 
IF     (N.EQ.NQ2)    GO    TO    603 
C 
C 

C  FCUR    POINT    LAGRANGIAN    INTERPOLATION 

C 

DJS=(T2*T3+T2*T4+T3*T4>*CJS(K,NQ2 1/6.0 

2  -IT3*T4+T4*T1+T1*T3)*CJS(K,NQP 1/2.0 

3  +<T4*T1+T4*T2+T1*T2)*CJS(K»NQ 1/2.0 

4  -(Tl*T3*T2*T3+T2*Ti)*CJS(K,NQM)/6.0 
0JPHI=<T2*T3  +  T2*T*+T3*T4)*CJPHKKfNQ2 1/6.0 

2  -<T3*T4+T**Tl+Tl*T3)*CJPHI(K,NQP)/2.0 

3  +<T4*Tl+T4*T2+Tl*T2)*CJPHI(K,NQ)/2.0 

4  -(Tl*T3*T2*T3+T2*Tl)*CJPHI(K,NQMl/6.0 
GO    TO      62  0 

602  CONTINUE 
C 

c 

C        TWO    POINT  LAGRANGIAN  INTERPOLATION 
C 

DJS   =  CJS(K,NQ  )-  CJS(K,NQM) 

DJPHI  =  CJPHI(K,NQ    1-CJPHKK,  NQM1 
GO    TO      620 

603  CONTINUE 
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C        THREE  POINT  LAGRANGIAN  INTERPOLATION 
C 

DJS=(T4+T3)*CJS(K,NQP)/2.0-(T2+T4l*CJS(K,NQ) 
2       +  <T2«-T3)*CJS(K,NQM)/2.0 

DJPHI=(T4+T3)*CJPHI(K,NQP)/2.0-(T2+T4)*CJPHI(K,NQ  ) 
2      •MT2+T3)*CJPHI(K,NQM)/2.0 

GO  TO   62  0 
601    CONTINUE 

IF  (N.EQ.NQP)  GC  TO  600 

IF  (N.EQ.NQ2)  GO  TO  604 
C 

C        THREE  POINT  LAGRANGIAN  INTERPOLATION 
C 

DJS=(T2+T3)*CJS(K,NQ2)/2.0-(T3+Tl)*CJS(K,NQP) 
2       +<T2+Tl»*CJS(K,NQ)/2.0 

DJPHI=(T2>T3I*CJPHI (K,NQ2)/2. 0-( T3  +  T1 »*C J  PHI ( K, NQP) 
2       +(T2+Tl)*CJPHI(K,NQ>/2.0 

GO  TO   62  0 
604    CONTINUE 
C 

C        TWO    POINT  LAGRANGIAN  INTERPOLATION 
C 

DJS   =  CJS(K,NQP)-  CJS(K,NQ) 
DJPHI=CJPHI (K,NQP  )-CJPHI(K,NQ) 
620    CONTINUE 
C 

C        CONTRIBUTION  OF  THE  K-TH  SPACE  SAMPLE 
C 

HSF(N)=HSF(N)*RO(K)*<SI(KI*DJS+DJPHI)*DS(K) 
C 

C        NEXT  SPACE  SAMPLE 
C 

600      CONTINUE 
650    CONTINUE 
C 

C        FINAL  VALUE  OF  THE  N-TH  TIME  SAMPLE  OF  THE  FIELD 
C 

HSF(N)=C^*HSF(N) 

X(N)=N*CP 

WRITE(6,9610)  N,HSF<N> 
9610  F0RMAT(5X,I5,5X,E12.5» 
C 

C  NEXT    TIME    SAMPLE 

C 

500         CONTINUE 
C 

C  END    OF    COMPUTATION    OF    THE    FIELD 

C  PLOT    SACK-SCATTERED    FIELD    VERSUS    TIME{ L I GHT-METERS) 

C 

WRITE(6»9630) 
9630       FORMAT(lHl) 

MODCUR=0 

CALL    PLGTP(X,HSF,ITS,MODCUR) 

IF     (INC.EQ-1)    GO    TO    112 

WRITE    (6,96241 
9624         FORMAT!//, 14X,'    NORMALIZED    BACK    SCATTERED    FIELD*,//, 

2  14X,'RAMP  RESPONSE', //,14X, 

3  'X-AXIS  TIME  (C*DT>» ,10X, 'Y-AXIS  FIELD  (H*C*R)») 
GO  TO  113 

112  CONTINUE 
WRITE  (6,9724) 

9724         F0RMAT(//,14X,'     NORMALIZED    BACK    SCATTERED    FIELD',//, 

2  14X, 'IMPULSE    RESPONSE', //,L4X, 

3  'X-AXIS    TIME     (C*0T1',10X, 'Y-AXIS    FIELD     (H*R)») 

113  CONTINUE 
RETURN 
END 
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Q  ************************ 

C  *    SUBROUTINE  SHAPE    * 

C  *                       * 

C  ************************ 

C 

C  PURPOSE: 

C  ******** 

c 

C  THIS  SUBROUTINE  COMPUTES  THE  SHAPE  GF  THE 

C  BODY  FROM  ITS  APPROXIMATE  COMPUTED  PHYSICAL-OPTICS 

C  RESPONSE. 

C  IT    FINDS    THE    CONTOUR    FUNCTION    RO(Z)     AND    DIVIDES    IT 

C  INTO    APPROXIMATELY    EQUAL    SEGMENTS    IDS). 

C 

C  USAGE: 

C  ****** 

C 

C  CALL   SHAPE  (  OS,  I  TER,  CP  , DSO  ,  IT  ) 

C  THE  DATA  NEEDED  AT  THE  INPUT  ARE  DSO, AND  HSF. 

C 

C  PARAMETER  DEFINITION: 

Q  ********************* 

C  ZMAX:  THE  MAXIMUM  LENGTH  OF  THE  BODY  ON  THE  Z-AXIS. 

C  (CORRESPONDING  TO  THE  POINT  WHERE  THE  RAMP  RESPONSE 

C  CHANGES  ITS  SIGN) 

C 

C  IT:  THE  TRANSIT  TIME  NEEDED  TO  OBTAIN  THE  REFLECTION 

C  FROM  THE  POINT  LOCATED  AT  ZMAX. (IT  IS  THE  NUM8ER 

C  OF  TIME  INCREMENTS). 

C 

C  DSO:  OPTIMAL  LENGTH  OF  A  SEGMENT (NEEDED  FOR  THE 

C  COMPUTED  SHAPE  SUBDIVISION)  .(METER) 

C 

C  ALL    OTHER    PARAMETERS     ALREADY    DEFINED    IN    MAIN    PROGRAM 

C 

C  REMARKS : 

Q  ******** 

c 

C  IN    THIS    PROGRAM    THE    RAMP    RESPONSE    USED    ASSUMES    A 

C  VERTICAL    H    FIELD    INCIDENT    (HORIZONTAL    POLARIZATION) 

C  THIS    GIVES    A    POSITIVE    INITIAL    RESPONSE. 

C  FOR    VERTICAL    POLARIZATION    THE     INITIAL    RESPONSE    IS 

C  NEGATIVE.     (THE   WHOLE    RESPONSE    CHANGES    ITS    SIGN) 

C  SINCE   THE    RADIUS    OF    THE    BODY    AT    EACH    POINT    IS 

C  PROPORTIONAL    TO   SQRT    OF    THE    INITIAL    RAMP    RESPONSE, 

C  IT    MUST    BE    ENSURED    THAT    THE    RAMP    USED   HAS    PROPER    SIGN. 

C 

C  METHOD  OF  SOLUTION: 

Q  ******************* 

c 

C  FIRST  STEP  IS  TO  FIND  ZMAX, WHICH 

C  CORRESPONDS  TO  THE  POINT  WHERE  THE  RESPONSE 

C  CHANGES  ITS  SIGN. 

C  THE  CONTOUR  FUNCTION  OF  THE  BODY  RO(Z)  IS 

C  PROPORTIONAL  TO  THE  SQRT  OF  THE  PHYSICAL-OPTICS 

C  RAMP  RESPONSE. 

C  THE  BOOY  IS  DIVIDED  INTO  APPROXIMATELY  EQUAL 

C  SEGMENTS, AND  THE  SPACE  SAMPLES  DELIMITING  THE 

C  SEGMENTS  ARE  FOUND  USING  A  »CUT  AND  TRY  METHOD' 

C  WHICH  CONVERGES  TO  THE  DESIRED  LENGTH  OF  SEGMENT. 

C  THE  RADIUS  AT  THE  SPACE  SAMPLE  POINTS  ARE 

C  INTERPOLATED   FROM  THE  KNOWN  SAMPLES  OF  THE  RAMP 

C  RESPONSE. 

C  THE  DATA  OBTAINED   Z ,RO,DS, NS, I S  USED  FOR  THE  NEXT 

C  COMPUTATION  OF  THE  FIELD. 
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c 
c 

SUBROUTIME       SHAPE    ( DS , ITER, CP ,DSO» IT) 
DIMENSION    DS(20)     ,RO(21),    Z( 21 ) , HSF( 60) 

COMMON    HSFt Z, RO,ITStNS,NSl 

COMMON    XMAX 
C 
C 

C        START  PROGRAM 
C 
C 

C        FIND  LENGTH  OF  BODY  (ZMAX  ) AND  ITS  CORRESPONDING 
C        TIME   OF  APPEARANCE  IN  THE  RAMP  RESPONSE  (IT) 
C 

DSP=DSO 

DO  710  N=l,  ITS 

IF(HSFdM)  .LT.O.O)  GO  TO  720 
710     CONTINUE 
720     CONTINUE 

IT  =N-1 

AD=1./(1.-HSF(N)/HSF(IT  )) 

ZMAX=(IT  -2+AD)*CP/2. 
C 
C 

C  COMPUTE    SHAPE    OF    BODY    UP    TO    ZMAX. 

C 

54  CONTINUE 
Z(l )=0.0 

RO    (1)=SQRT(2*HSF(2) ) 
DO    50    K-1,19 
K1  =  K+1 
U=DSP/2.0 

55  CONTINUE 
NN=1 

C 

C  ADJUST    Z    TO    OBTAIN   PROPER    DS 

C 

Z(K1)=Z(K)+U 
53  CONTINUE 

IF(Z(K1).GT.ZMAX)G0    TO    56 

GO    TO    57 

56  CONTINUE 
Z(K1)=ZMAX 

RO(KU  =  0.0 
GO    TO    51 

57  CONTINUE 
C 

C  COMPUTE    CORRESPONDING    RADIUS    ( RO >    OF    THE    BODY 

C 

G=2*Z(K1)  /CP+2. 

M=IFIX(G) 

M1=M+1 

RQ2  =  (HSF(M1)-HSF(M)  )*( G-M)+HSF( M ) 

R0(K1)=SQRT(R02*2.) 
C 

C  COMPUTE    LENGTH    OF    SEGMENT    DS 

C 
51  CONTINUE 

DS     (K)  =  SQRT((Z(K1)-Z(K)  )**2  +  ( R0( K 1  )-R0( K)  )**2  ) 

A=DS(K)/DSP 
C 

C        ALL  OBTAINED  SEGMENTS  MUST  BE  WITHIN  1/100  OF  DSJ 
C 

IF  (A.GT. I. Oil  GO  TO  58 

IF( Z(K1 ).EQ.ZMAX)GO  TO  60 

IF  (A.LT.0.99)  GO  TO  58 

GO  TO  50 
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c 

C  ADJUST    Z    AND    GD    BACK    TO    CHECK    FOR    DS 

C 

58  CONTINUE 

NN=NN+1 

IF     (A.GT. l.OJ    GO    TO       90 

Z{K1)=Z(K1)+U/NN 

GO    TO    53 
90  CONTINUE 

ZiKl)=Z(Kl)-U/NN 
IF    (Z(Kl).GT.Z(K))     GO    TO    53 

Z(K1)=Z(K1)+U/NM 
NN=NN+1 

GO    TO    90 


50 

C 

C 

C 

60 

CONTINUE 

NEXT    SEGMENT 

CONTINUE 

NS1=K1 

C 

c 
c 

NS=NS1-1 

ADJUST    FOR    LAST    SEGMENT 

U2=DS(NS) /DSP 

IFCU2.LT. 0.10)    GO    TO    62 

IF<U2.LT.0.95)    GO    TO   61 

GO    TO    70 

61 

CONTINUE 

U3=(Z<NSU-Z(NS))/(NS-l) 

NS1=NS1-1 

NS=NS1-1 

DO    80    K=i,NS 

K1=K+1 

ZIK1)=Z(K1)+K*U3 

Z(NS1)=ZMAX 

G=2*Z(K1) /CP+2. 

M=IFIX(G) 

M1=M+1 

R02=(HSF(M1  )-HSF(M)  )*(  G-M  )+HSF(  M) 

R0(K1)=SQRT(R02*2.) 

RO(NS1)=0.0 

DS    (K)=SQRT<< Z(K1 )-Z(K)  )**2  +  ( R0( K 1 )-R0< K ) )**2  ) 

30 

CONTINUE 

GO    TO   70 

62 

CONTINUE 

Z(NS)=ZMAX 

RO(NS)=0.0 

NS2=NS-1 

DS(NS2)=SQRT( (ZMAX-Z(NS2) ) **2+R0( NS2)**2) 

NS1=NS1-1 

NS    =NS1-1 

70 

c 

c 
c 
c 
c 

CONTINUE 

SHAPE    OF    BODY   COMPLETED    AND    SUBDIVIDED 

WRITE    NUMBER   OF    SEGMENTS    OBTAINED 

WRITE    (6,9601  )    ITER,NS 
960L       F0RMAT(1H1,5X,»NEW    NUMBER    OF    SEGMENTS    AFTER     ITERATION* 
2       tI4,«     IS» ,I^,///,10X,,Z,,10X,«RO« ,10X,«DS» ,/////) 

DS    (NS1)=99.9 

DO    40    N=l,NSl 

WRITE(6,9602)    Z<  N  )  ,  ROC  N)  ,DS    (N) 
9602  FORMAT(5X,3F10.5I 

40  CONTINUE 

C 
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c 

C  PLOT    THE    CONTOUR    FUNCTION    OF    THE    BODY 

C 

NS2=NS1+1 

Z<NS2)  =    XMAX 

RO    (NS2I=    Z(NS2) 

WRITE(6,9603) 

9603  FORMAT* 1H1) 
MODCUR=0 

CALL    PLOTP( Z,R0,NS2,M0DCUR) 
WRITE    (6,9604)    ITER 

9604  FORMATC// ,20X, 'SHAPE    OF    BODY    AFTER    ITERATION    NUMBER    », 
2    I4,//,20X, »X-AXIS    Z      Y-AXIS    RADIUS    (R0)») 

C 

C        END  OF  SUBROUTINE  SHAPE 

C 

RETURN 

END 
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c 
c 

c 
c 

C  EXAMPLE    OF    DATA    CARDS 

C 

C 

c 

0.0  0.0 

0.02185  0.20791 

0.08645  0.40674 

0.19098  0.58779 

0.33087  0.74314 

0.50  0.86603 

0.69098  0.95106 

0.89547  0.99452 

1.10453  0.99452 

1.30902  0.95106 

1.50  0.86603 

1.66913  0.74314 

1.80902  0.58779 

1.91355  0.40674 

1.97815  0.20791 

2.0  0.0 
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APPENDIX  F 
RELATIONSHIP  BETWEEN  RAMP  RESPONSE  AND  SHAPE 

The  validity  of  (5.2)  will  be  demonstrated  analytically 
here  in  the  case  of  a  metallic  body  of  revolution.   We  first 
express  the  axially  directed  incident  ramp  field  as: 


H^(z,t)   =   {  (F.l) 


t  - z/C     t  >  z/C 
t  <  z/C 
The  physical  optics  current  induced  on  the  body  is  given  by: 


J(r,t)   =   2a  xS(r,t)  (F.2) 

n    R 


The  backscattered  far- zone  field  is  given  by  equation  (3.8) 
and  its  scalar  equivalent  for  the  body  of  revolution  is 
expressed  by  (3.16).   Substitution  of  the  current  in  (3.16) 
gives  the  physical  optics  ramp  response  as  follows: 


H|(t)   =   "4^0   /   |^[-4H^(z'/T)sin  a'lp'ds'     (F.3) 

The    time   derivative    of   the   incident   field  expressed  by    (F.l) 
is   unity: 

TF=     TF    =     1         t<z^c  (F-4) 


Substitution  of  (F.4)  into  (F.3)  gives: 
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s         1    .s 

HR(t)   =  FT  J   sin  a'  p'ds'  (F>5) 

o   0 


From  the  geometry  of  the  body  of  revolution  we  note  that: 

dp'   =   sin  a1  ds' 

Substitution  in  (F.5)  gives  the  following  expression: 

p(z) 


(t)   =  -i-  ./   "  p'dp'   =   2^c  p2(z)      (F-6) 


o   0 


Since  the  area  function  is  given  by: 


2 

A(z)    =    TT  p   (z) 


the  physical  optics  ramp  response  (F.6)  will  be  related  to 
A ( z )  by 


HR(t»   "   2^-CA(Z>  (F-7» 
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